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ABSTRACT 





This-thesis examines the-system implementation considerations of a resonance- 
based. radar target classification system. The basis of the system is the aspect and 
excitation independent property of electromagnetic scattering from a conducting 
body. Such a system consists of two components: pole extraction and annihilation 
filtering. The algorithms investigated here for these purposes are the Cadzow- 
Solomon pole extraction-algorithm-and the K-Pulse annihilation filter. Additionally, 
an aspect-dependent annihilation filter based on an inverse ARMA model is 
introduced. The procedures are applied to noise polluted synthetic data, as well as 


scattering data collected for a thin-wire and silver coated 1/72 scale mode] aircraft. 
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J. INTRODUCTION 


Modern antiair-warfare has, with the advent of radar in World War II, become 
a process of detecting aircraft and: identifying them as friend or foe at sufficient 
ranges to allow engagement with appropriate weapon systems. Radar development 
has historically emphasized the first.element of this process: detection. As modern 
air warfare has become more complex, dynamic, and lethal, identification has become 
as critical as detection. 

In the past, identification has been facilitated-by active measures on the part of 
friendly aircraft to distinguish them from enemy, i.e. IFF transponders, specific flight 
routes, etc. All these are subject to exploitation by an adversary. This has given rise 
to development of non-cooperative target recognition (NCTR), the ability to identify 
an enemy aircraft despite its efforts to deceive. 

One approach to NCTR, as discussed in this thesis, is identification based on 
the natural complex resonances of the aircraft as a scatterer of incident 
electromagnetic (radar) energy. In 1971 Baum showed that an_ incident 
electromagnetic wave creates a current response on a body with certain complex 
natural resonances dependent on its physical configuration [Ref. 1]. Moffat 
and Mains proposed that a target scatters electromagnetic energy with resonances 


that are independent of incident polarization or angle [Ref. 2]. Morgan 


confirmed this through a rigorous theoretical description of induced currents and the 

















resultant scattered: field: [Ref..3]. The problem, then, is to determine these 
complex poles for a target of interest and ‘subsequently use -them:.to identify this 
target. 

The first step in using complex natural resonances for target identification: is 
extraction of these poles from scattered responses of actual targets (either modeled 
-or full size) [Ref..4]. These poles are then used to create annihilation filters 


for each potentially observed target [Ref. 5]. A given filter will, in principle, 


produce a near null late-time response to its matched target’s scattered signal,. 


independent-of aspect or polarization. 

The objective of this thesis is to.examine these processes in order to define the 
overall system characteristics -and operational considerations of a radar target 
identification system based target natural .resonsnces. This analysis will use the 
application of signal -processing algorithms to -both synthetic and actual 
electromagnetic scattering data. This-thesishas been.conducted in-conjunction with 
research in pole extraction by Murphy [Ref. 6], whereby the data generated 


in that work was incorporated in the filters discussed in -this effort. 


A. BACKGROUND 

The ability to extract poles and ‘use them to identify a target is dependent-on 
the nature of electromagnetic scattering from a conducting body. Consider a 
generalized electromagnetic pulse with electric and magnetic field components E'(r,) 


and H'(r,t) respectively, travelling in a direction. p and incident upon a finite-sized, 








perfectly conducting body in.free-space. Currents-are induced on-the body which are . 


described by the magnetic field integral equation: [Ref. 7], 


I(rt) = WxH'(,p + f ie Reese's Phas (1) 
PY 


where n is the outward unit-normal vector, J is the surface current density, and K is 
a dyadic Green’s function kernel. Here the surface integral is a principle-value type, 
excluding the point r = r’. The cross product term represents-the physical optics 
contribution to the induced currents and the surface integral represents the feedback 
of:currents on each point on the object-due to the currents-at.all other points. This 
feedback mechanism forms the physical basis for natural resonances. 

After the incident pulse-leaves the target body, H’ in Equation (1) goes to zero; 
the source-free current modes that remain-"ring" at the natural_resonance frequencies 
of the object. These currents are described by the expression J Ar)exp(s,,t); where 
S, = 6, + jw,, are the frequencies which are dependent only on the physical 
properties of the scatterer. These complex frequencies are described in the Laplace 
domain by their damping constant, ¢, in Nepers/sec, and their frequency, o, in 
radians/sec. An infinite number of these resonant modes are inherent in any 


candidate scattering object, but only those in. the bandwidth of the incident pulse are 


significantly excited. [Ref. 7] 




















The induced currents on the scattering ‘body, J, give rise to -a far-field: 
distribution-in the backscatter-direction, -p, (the direction of-interest-for monostatic 


radar applications) described ‘by 


Srp) = 2 26 perro has! ® 
H-p) = 7 Sf, Pel as’ 


This equation is applicable during the entire time-the incident pulse-is traversing the 
target.in the direction p: as well as.after it has gone beyond it. When: the specific 


current distribution described in Equation (1) is:substituted here, the result is 


HS(-1p,) = u(t-r/c) | H,(-9.) + YD H-rBdexp(s;p (3) 


a~co 
nvO 


where H,, is a physical optics term dependent on the incident direction, p, relative - 
to the object shape, and the summation term is due to the source-free currents 
‘behind the advancing wavefront. The field-due to these currents. resonates. at the 
same frequencies as the currents themselves. However, during the time the pulse still 
resides somewhere on ‘the target, the coefficients H,, are varying. After the pulse 
leaves, the aspect dependent physical optics term vanishes and the coefficients of the 
scattered resonant frequency modes become constant. This instant defines the 
transition from “early-time" to "late-time" in the scattered field. For backscattering 


this transition occurs after a time given by T, = 2D/c + 7, as measured from when 








the pulse first impinges on the target, where D i: the length of the object in the 
‘propagation direction, j; c, the speed of light; and 7, the incident pulse width. [Ref. 


q 


B. PROBLEM 

Application of these characteristics to target identification involves two separate 
steps: pole extraction and annihilation filtering. The transfer function of an 
electromagnetic target contains a complicated set of poles that can be analytically 
predicted only for the simplest objects, e.g. a thin-wire or a sphere [Ref. 8]. 
Therefore, it will be necessary. to extract the poles from complex targets, such as 
aircraft, through signal processing of their scattered response [Ref.4]. A number of 
algorithms have served as candidates for this process [Ref. 4], including the 
Kumaresan-Tufts and. Cadzow-Solomon algorithms compared by Larison in 
[Ref. 9]. These will be discussed further in Chapter IIL. 

After the poles are determined for the targets of interest, they will be 
incorporated in a bank of annihilation filters; each filter corresponding to a target. 
Scattering responses from unknown targets (radar returns in actual application) are 
passed to all the filters. The annihilation filter with the least output should correlate 
the unknown target to a previously cataloged response, thus identifying that target. 
This is the basis of the Resonance Annihilation Filter (RAF) concept discussed in 


[Ref. 7] and [Ref. 10]. Further analysis of this process and its implications- 


for actual implementation are discussed in Chapter IV. 








€. HISTORY 

As stated, the.initial workin the area-of natural-resonance response was derived 
from-Baum’s investigation-of currents generated as a:result of electromagnetic:pulse 
(EMP) [Ref. 1]. ‘Fesche developed the-concept further into the specifics of the 
Singularity Expansion Method:(SEM).applied to-generating th time response of a 
thin wire scatterer analytically by examining its poles and residues [Ref. 11}. 
Morgan-points out in [Ref. 3] and (Ref. 7] that the SEM description-of the scattered 
signal is-applicable only in late-time, the time after the incident pulse leaves the 
object. As a result of this limitation, earlier -schemes for pole extraction based on 
Prony’s method [Ref.12] were limited by noise and clutter [Ref. 13]. 

Since that time, effort has been directed by many researchers into improved 


pole extraction algorithms. Chong incorporai2d-the early-time in using a direct non- 


linear least-squares minimization [Ref. 14], >dut this approach is. 


computationally intensive and requires much a priori knowledge of the signal [Ref. 
4]. The Kumaresan-Tufts method of pole extraction-is limited in applicability to the 
late-time only, but relieves some of the burden of # priori knowledge of the natine 
of the scattered response [Ref 15] and improves the noise handling 
capability of Prony’s method. The Cadzow-Soiot... 2 technique, a different extension 
of Prony’s method, makes use of both the e?;ly and late-time, allowing it to use as 


much of the available energy as desired in the return signal. This algorithm, 


introduced in [Ref. 16] and applied to pole extraction in [Ref. 4], [Ref. 9] 


and [Ref. 6], does require input of the transmitted signal along with the return 























scattering signal, but holds the most promise of the considered extraction fnethods 
in the presence of noise. 

The annihilation filtering concept, #3 described above, has been attempted 
primarily with two basic algorithms. ‘vy: » alse algorithm, based on-a difference 
equation representation of the filter, kins "s <n developed by Morgan [Ref. 7} and 
Dunavin [Ref.17] from the  oripy.al work -of Kennaugh ([Ref. 18}. 
The E-Puise technique, developed by «+ i [Ref. 19], uses. integral operators 
to cancel the resonances.in late-time. Simun compared these algorithms in [Ref. 1¢: 
and found similar results. This thesis will further investigate the K-Pulse annihilation 
filter, and will extend it to.a recursive versicn which works on early-time energy as 


well. 














TI. EXPERIMENTAL. SIGNALS: 

Both synthetic-and:actual signals were used-to analyze:the per;crmance of:pole 
extraction and annihilation filtering. Two types of synthetic:signals were used. First,. 
synthetic signals with predetermined complex poles were used to val! 2ate.extraction- 
an. annihilation software end to gain-insight into the mechanisms by which these 
algorithms worked. Secondly, poles were extracted from synt. ctically generated:thin- 
wire scattering data generated by Morgan’s thin-wire integral equation program: 
Ref. 20], and ieee in annihilation. filters. For this set of data, the poles. 
weré not -known in advance. After the software was.tested on synthetic data, it was 
applied to actual scattering data, gathered by Walsh [Ref. 21]. This data was 
measured for a thin. wire and several scale model aircraft using the NPS Transient: 
Electromagnetic Scattering Luooratory:(TESL). In the development.of the first set 
-of synthetic signals, as many parameters as-possible were selected to. izatch those 


anticipated in the TESL measured data. 


A. ARMA SYNTHETIC SIGNALS 


]. 'Scatterer Model 
As discussed in the previous chapter, the scaitered field consists of a 
physical optics term due to the driving incident wave and a feedback portion 
representing the complex natural resonances. These characteristics suggest an. 


autoregressive moving average (ARMA) or recursive structure for modeling this 











ee 








signal. Such a structure was realized with a recursive difference equation signal 
generator, 


k=N k=L- - 
yn) = Vapn-b + Yebx(n-h) @) 
k=0. 


kel 


where x(1):is the digital version of the incident pulse, y(7) is the scatterer response 
and-a, and-b, are coefficients determined hy the-complex natural resonances and 
initial conditions, respectively [Ref.22:pp. 81-84]. After selection of a 
sampling frequency-for the syntsetic signals, the coefficients were determined through 
conversion of a simple set of .-plane poles to the complex z-plane, followed -by 
multiplication of rational polynomial terms representing. each-pole pair to arrive at 
an overall ratio of polynomials: the z-transform version of the scatterer’s transfer 
function. The coefficients from the numerator and denominator appear in Equation 
(4) as b, and a,, respectively. Additional inputs were a complex coefficient of the z- 
plane poles to allow.scaling or phase differences betweer: poles. Program listings for 
the coefficient generator and_the recursive signal generator appear in Appendices A 
and B. | 

The poles were selected in the s-plane to coincide generaily in-bandwidth 
with those obtainable by actual scattering measurements. As such, three-sets of ten 
pole pairs each were selected, lying at integer values between 1 - 10 GHz. The three 


sets represent aifferent levels of damping (Q factor) and.are related in that they all 








have:the same-modulus when plotted as-vectors in the s-plane. For each set of poles, 


the frequency-and:damping coefficient, w, and a, are. related-by a:damping factor k,, 


Go 


= == [ink] (5): 
where 

k, = 0.5.(low Q/highly damped) 

k, = 0.7-(medium Q/moderately damped). 

k; = 0.8:(high Qflightly damped) [Ref. 10]. 
These poles are plotted in the-s. and z-planes respectively ‘in Figures 1.and:2 and are 
listed in Tables 1, 2 and 3. The units of f,, 9, and w, are GHz, GNpj/sec and 


GRad/sec, respectively. The resulting ARMA coefficients, a, and-b,,.are in Table 4. 


TABLE 1. LOW Q SYNTHETIC POLES 


POLES RESIDUES 

Sp op G, 4, p 
1 0.6892 6.2474 1 0 
2 1.3784 12.4948 i 6 
3 -2.0676 18.7422 1 0 
4 2.7568 24.9895 1 0 
5 -3.4460 31.2369 1 0 
6 4.1352 37.4843 1 0 
7 4.8244 43.7317 1 0- 
8 5.5136 49.9791 1 0 ? 
9 6.2028 56.2264 1 0 
10 68919 62.4738 1 0 
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Figure 2. Synthetic poles, z-plane 
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TABLE. 2. MEDIUM Q SYNTHETIC POLES 











POLES RESIDUES 
P op p A, ®p 
1 0.3562 -6.2752 1. 0 
2 0.7124 12.5504 1 0 
3 “1.0687 18.8256 1 0 
4 -1.4249- 25.1007 1 0 
5 “1.7811 31.3759 1 0 
6 -2.1373° 37.6511 1 0- 
7 2.4935 43.9263 1 0- 
8 -28498 50.2015 1 0- 
° -3,2060 56.4767 1 0 
10 3.5622 62.1518 1 0 

7 TABLE 3. HIGH Q SYNTHETIC POLES 

- POLES RESIDUES 
Sp op Op A, Vp 
1 0.1625 6.2832 1 0 
2 0.3250 12.5664 1 0 
3 0.4876 18.8496 I 0 
4 0.6501 25.1327 1 0 
5 0.8126 31.4159 1 0 
6 0.9751 37.6991 1 0 
7 “1.1376 43.9823 I 0 
8 +1.3002 50.2655 1 0 
9 -1.4626 56.5487 1 0 
10 -1.6252 62.8319 I 0 











The sampling frequency of 25.6 GHz, 512 samples in 20 nanoseconds, was. 


‘used in the creation of the synthetic signals. This is- the minimum "power-of-two" 





TABLE 4. ARMA COEFFICIENTS :-FOR: SYNTHETIC SCATTERING 





Low Q Medium Q High Q- 
k &, 4%, « b, “; b, 
0 10 20.0 1.0 20.0 1.0 20.0 
1 “3.77 “71.59 “3.54 67.32 -3.46 65.65 
2 8.47 152.61 8.08 145.36 7.96 143.26 
3 -14.27 -2A2.52 -14.09 -239.57 “14.14 -240.43 
4 20.02 320.34 21.00 335.97 21.66 346.64 
5 -24.51 -367.64 -27.61 414.21 -29.40 441.01 
6 27.07 379.04 33.15 464.08 36.58 $12.05 
7 -27.40 -356.18 -36.69 476.94 42.03 -546.42 
8 25.80 309.55 38,02 456.23 45.32 543,92 
9 -22.70 -249.71 -36.96 -406.56 ~45.90 -504.90 
10 18.84 188.36 34.00 339.98 44,02 440,24 
11 -14.72 -132.44 -29.50 265.52 -39,85 -358.65 
12 10.89 87.13 24.28 194,28 34,22 273.76 
13 -7.58 -53.04 -18.80 -131.63 -27.64 -193:49 
14 4.99 29.92 13.75 82.51 21.10 126.40 
15 -3.04 -15,22 9.33 46.66 “14.89 -74.45 
16 1.74 6.98 5.90 23.60 9.77 39.09 
17 -0.89 -2.67 -3.33 9.99 5.72 17.17 
18 0.42 0.84 1.69 3.39 2.99 5.98 
19 0.15 0.15 0.67 0.67 “1.21 “1.22 
20 0.05 0.00 0.22 0.00 0.38 0.00 


numvuer of samples allowable in a 20 nanosecond collect time by the Nyquist 


requirement applied tothe frequencies of interest, one.to ten GHz. The requirement 
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to limit the sampling bandwidth as much-as-possible is derived from the use of the 
annihilation-filter, discussea- Chapter IV. 
2. Incident Pulse 

After determination of the coefficients of the synthetic scatterer. transfer 
function, these were used in-the difference equation, Equation (4), above. The input 
waveform, x(n), was chosen to be a double Gaussian, in order to match that-used in 
the actual scatterer data collection. The double Gaussian waveshape consists of a 
broad, slow Gaussian pulse-subtracted- froma narrow, rapid Gaussian pulse. The 
narrow pulse determines the high frequency roll-off; the broad pulse determines the 
low frequency roll-off. The incident pulse-used here has ten percent narrow and wide 
pulse widths of 0.15 and 0.3 nanoseconds, respectively, and-is pictured in Figure 3. 
The frequency domain representation is in Figure 4. The final results of the double 
Gaussian input to the ARMA model is shown for all three pole sets (high Q, sbediaii 
Q, and low Q) in Figures. 5, 6, and 7 and Figure 8 contains the spectrum of the 


medium Q target. 
3. Noise Modeling 
The synthetic signals were polluted- by simulated Gaussian noise at various 
signal-to-noise ratio levels: 20 dB, 15 dB, 10 dB and 7 dB. The noise was introduced 
by a noise pollution routine, Appendix C, which uses a random number generation 


algorithm based on subroutines of the Fortran Scientific Library. This routine 


generates random numbers under a uniform distribution. The Central Limit 
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Figure 3. Incident double Gaussian pulse 
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Figure 4. Frequency spectrum of double Gaussian 
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Figure 5. High Q synthetic scattering 
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Figure 6. Medium Q synthetic scattering 
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Figure 7. Low Q synthetic scattering 
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Figure 8. Frequency spectrum of medium Q synthetic target 
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Theorem was then invokeu by calling this routine 1600 times for each-additive noise 
point to ensure the-resultant noise was Gaussian distributed. {Ref. 6]: Figures 9 and 
10 show thée-high Q waveform of Figure 5: with the:addition of-hoise of 20 dB and'7 


dB, respectively. 


B. THIN-WIRE INTEGRAL-EQUATION. SIGNALS. 

The concept of extracting-unknown poles and then incorporating them into 
filters was validated on the signals generated:by the Time-Domain Integral Equation 
{(TDIz) Thin-Wire Scattering Program, due to Morgan [Ref. 20]. This:program uses 
a record length of 960 samples in a 20 nanosecond collect time. The scattering 
simulated excitation by the same double Gaussi?n-pulse-with approximately a 1-10 
GHz bandwidth. The program-generates wire.scattering data for any combination 
of incident and reflection:angles. Backscatter data, generated at angles of 175, 150, 
110, and 90 degrees measured relative to the long axis of the wire, were used to 
effect a monostatic situation. This data was similarly noise polluted at the same 
signal-to-noise ratios. A-representative case-for the noiseless TDIE thin-wire data-is 


in Figure 11 and its spectrum is in Figure 12. 


C. LABORATORY: SIGNALS 

Actual’scattering signals used for. pole extraction and annihilation filtering were 
obtained from the NPS Transient Electromagnetic Scattering Laboratory. As stated, 
the parameters determined by the instrumentation of the laboratory were used. 


throughout the experiments for both synthetic_and-actual signals as much as possible. 
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Figure 9. High Q synthetic scattering with 20-dB SNR 
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Figure 10. High-Q synthetic scattering with 7 dB SNR 
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Figure 11. TDIE thin-wire scattering with 30° incidence, 150° reflection 
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Figure 12. Spectrum of TDIE thin-wire 367150° data 
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This included the scattering collect time of 20 nanoseconds. and’ the -excitation 
bandwidth-of 1-10: GHz due-to the:incident pulse-:double Gaussian -waveshape. The 
sampling frequency.serves as the notable exception to this matching. The TESL 
samples 1024 times in the 20:nanosecond collect time for-a sampling frequency of — 
51.2:GHz. This twice oversampling, while allowing more flexibility in the-use-of the 
TESE, dees as-a limiting factor for subsequent annihilation filtering. As a-result, the 
collected scattering data from the Jab was collapsed to the 25.6 GHz sampling 
‘frequency, mentioned above, by averaging each-pair of points in the record before 
annihilation filtering. The TESL provides a. nominal signal-to-noise ratio of 


approximately 26 dB [Ref. 21]. 


1. Deconvolution Postprocessing 

The double Gaussian waveshape, as described above, is not actually 
transmitted in the TESL. Instead, a noise-like -but ‘highly repetitive signal, with 
approximately the same bandwidth is transmitted and the response of the target to 
the specific double Gaussian waveshape is estimated through deconvolution post- 
processing. This allows the target to be subjected to an analytically precise-signal like 
‘the double Gaussian without requiring a tremendously sophisticated transmitter 
apparatus. Additionally, the free-space response of the target is more nearly 
obtained by removing the effects of the chamber and antennas. The details of this 


process are in (Ref. 21]. Targets used in the scattering experiments included a 0.1 


meter thin-wire target and silver coated 1/72 scale model aircraft. 

















2. Measured Thin Wire Target 

The transition to real-world scattering data-was.first made with scattering 
data coilected in the TESL on a 0.1-meter long, 2.36 millimeter diameter thin metal 
wire [Ref. 21]. This wire data was collected monostatically-at 30 degrees, 45 degrees 
and broadside to the wire. The time and frequency domain data for the 30 degree 
case are in Figures 13-and 14. Notice that the FFT of this target’s scattering is much 
more complicated than for the synthetic TDIE thin wire. The spectrum in Figure 14 
displays.a significant resonance at approximately 3 GHz, with less excited poles at 1.5 
and 4 GHz. As the apparent length of the wire, D, at this angle is 0.05 meters, so 
these resonances correspond roughly to modes at 4 = 2D, 3D, and 1.5D, respectively. 
This demonstrates the complicated interaction between the excitation and resonances. 


of the target. 


3. Model Aircraft Targets 
The bandwidth of the TESL, 1-12 GHz, is matched to the scaling factor 
of the model aircraft targets. Table 5 contains the significant dimensions, both actual] 
and scaled, of the aircraft used to collect scattering data. Comparison of the scale 
model lengths to the wavelengths of the TESL excitation, 2.5 to 30 cm, indicate that 
large fractions of wavelengths can be contained in the dimensions of the aircraft 


models. This allows for the lowest order, i.e. the strongest available, modes to be 


excited in the model structure. 
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Figure 13. Measured thin-wire Scattering at 30°/150° 
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Figure 14. Spectrum of 30° measured. thin 








TABLE 5. DIMENSIONS OF AIRCRAFT TARGETS [Ref. 23,.24] 


Target number 7 1 20 4 ss 

Overall length ; actual (m) 19.43 15.03 1707 12.20 15.88. 
model (cm) 26.97 20.88 23.71 16.94 22.06 
Wingspan actual (m) 13.05 9.45 11.43 10.96 13.95 
mode? (cm) 18.13 13.13 15.88 15.22 19,38 


The aircraft.data was collected by Walsh at six different aspects: nose-on, 
tail-on, broadside, 30 degrees-off the nose, and from above with both the fuselage 
and wingspan parallel to the incident electric field [Ref. 21]. Target No. 2'scattering 


data at thirty degrees aspect angle is in Figure 15 with-its spectrum, Figure 16. 


D. SUMMARY 

Each of these target sets; the ARMA synthetic scatterer, the TDIE thin-wire 
scatterer, the actual thin wire, and the model aircraft; each represent an increasingly 
more complicated and challenging case for the pole extraction and annihilation 
filtering algorithms. This is especially evident when reviewing their various frequency 


response plots. 
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Figure 15. Scattering from aircraft target No. 2, 30° aspect 
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Figure 16. Spectrum of target No. 2, 30° aspect 
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II. POLE:EXTRACTION 


The pole extraction algorithm considered in-this thesis is due to-Cadzow. and: 


Solomon [Ref. 16]. Other. algorithms have been used and are briefly discussed in. 


‘Chapter I. More detailed analysis of the other algorithms can be found in [Ref. 4]. 


A. CADZOW-SOLOMON ALGORITHM 

The scattered signal structure discussed.in Chapter I and dispiayed in Chapter 
II consists of an early-time-and-late-time separated: by the point where the incident 
signal ceases to excite the scatterer. In early-time, the resonances are just being 
initialized, so their amplitudes are small compared-to the driver. However, the end 
of early-time marks the point where they have their largest instantaneous power. The 
Cadzow-Solomon algorithm-takes advantage of this fact by attempting pole extraction 


from a data record containing both early and late-time data segments. 


1. Applicability 
The Cadzow-Solomon method seeks to identify the unknown coefficients 
a, and b, for an autoregressive moving average (ARMA) process given by Equation 


(4), which is repeated here, 


k=N k=L 


yn) = Dayn-b + Ybx(n-h, (6) 
ks0 


k=} 
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by analyzing the system input, x(n) and output, y(z) [Ref. 16]. Identification-of these 
coefficients will allow for extraction of the poles and zeros of:the transfer function 


which: has been acted on: by x(11) to ‘produce y(7), specifically, 
(7) 


2. Methodology 

The algorithm: does-not require specific knowledge-of the system order, NV 
and L in Equation (7), as it uses an "over-estimation" of this-order to account for 
noise and extract the poles. Norton.applied a technique used in the Kumaresan-Tufts 
algorithm, whereby the relationship in Equation (4) is written in reverse time order. 
This forces the desired poles to be-extracted outside the unit circle on the z-plane, 
while the extraneous poles, introduced by "over-estimation" to allow for noise, are 
constrained to lie inside the unit-circle [Ref. 4]. This reverse time, non-causal form 


of Equation (4) is 


k=N k=L 
y(n) = Say(Nen-k+1) + Yo ba(n-h. (8) 
k=} k=0 


The pole locations are then reflected through the unit circle resuiting in the signal 
poles within, noise poles outside. This final presentation of the poles describes the 
actual scatterer: a stable, causal system. 


For M segments of input-output data related by Equation (8), where 


















Pur “ Yew 7 % Yt 
i me 3 , 4 (9) 
Voom" Voenem-r * mer | Vesa | 
where both parts of the y|x data matrix are Toeplitz matrices in which each row é 


introduces one new value and shifts all others one placeto the right. Rewriting this 


equation in more compact notation yields 


Ds = y > (10) 


Extraction is accomplished by solving this matrix-vector equation for the ARMA 


parameter vector, 


Ci, = Dt y. (11) 
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where C,,* is the-minimum Euclidean norm-solution vector for C,, and D,,* is the 
‘psuedoinverse of the data matrix Dye The solution vector, C,,*, is the minimum 
“ Jength least-squares solution to the matrix-vector-equation. Refs. 4-and. 9] 
3. Singular Value Decomposition 
The above solution relies on singular value decomposition (SVD) of the 
‘data matrix, D,,, to derive its psuedoinverse, D;,". This process separates D,, into the 


product of three matrices, 


D. = Uzy" (12). 


-where U and F are both orthogonal matrices. (that is, U'U=) and where " denotes 
the transpose of V. The columns of the matrices U (of size M by M) and V (N by N) 
are composed of the eigenvectors of D,,D,,’ and D,,'D,, respectively. The singular 
values of the data matrix are the nonzero components of the M by N diagonal matrix, 
%. These singular values are the nonzero eigenvalues of both D,D,,” and D,,"Dy. 


Rearranging these components of the SVD of D,, produces its psuedoinverse, 


= anyT 13 
Dy, = VE"U". (13) 


Here, the diagonal components of &* are the inverse of the original singular values 


[Ref. 9]. 
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The original singular values must be compensated:for noise if it is:present 
in the input/output data which comprises. the data matrix. This compensation 
assumes the input and output noise variances to be equal. More detail on bias 


compensation is provided in [Refs. 4 and 9]. 


B. PERFORMANCE 

The Cadzow-Solomon algorithm:-has demonstrated the ability to-extract poles 
from scattering data. Larison [Ref. 9] applied the algorithm:to synthetic target data, 
both measured and analytically predicted thin-wire scattering-data, and aircraft model 
data. Murphy [Ref. 6] refined the process by further-application to synthetic data 
sets as well as measured thin-wire and-aircraft-data as presented in Chapter II of this 


thesis. 


1. Noise Tolerance 

Murphy tested the synthetic scattering at various signals-to-noise ratios 
down to seven dB. Figure 17 shows that the algorithm-remained relatively robust, as 
determined by observing the tightness of the clustering of the extracted poles, down 
to 20 dB [Ref. 6]. Below this level the poles would not cluster significantly to lend 
a high confidence to their location. Figure 17 presents fourteen extraction trials of 
medium Q data, each containing a different.20 dB noise sequence . Notice that the 
higher order poles (near 7/2) are not extracted. Fewer poles are available with 


successively greater noise levels [Ref, 6]. 
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Figure 17. Extracted poles, medium Q synthetic data, 20 dB SNR [Ref. 6] 











2. Accuracy Versus Precision 

The adequacy of the extracted poles can be couched in.terms of accuracy, 
the:distance ofthe extracted poles to:their true location; and_precision, the tightness 
of a group of poles, regardless of its location. Figure 17 demonstrates some 
degradation in accuracy as the extracted poles are slightly. biased toward the unit 
circle [Ref. 6]. 

The-noiseless TDIE data yielded excellent results in-producing very tightly 
grouped poles. This was to be expected because-of the-noise-free environment and 
the high-quality resonator represented by the modeled thin wire. 

For actual scattering data, accuracy can only be determined through 
evaluating the target discrimination effectiveness of the annihilation filters constructed 
from the extracted poles. Figure 18 shows a scatter plot of extracted poles from the 
measured thin-wire scatterer. The poles are represented by a scatter because their 
position can be perturbed by noise and is highly sensitive to the parameters used in 
the extraction algorithm. As such, the extracted poles from actual targets will 
represent some distribution about their true (unknown) location. Murphy used an 
averaging scheme on the extracted poles which lie in tight groups, after discarding 
outlying positions, for a given angle of incidence. These averages were then averaged 
over the available angles of incidence in an attempt to ameliorate any possible bias 


due to the parameters and to compensate for the changing relative excitation of one 


mode over another with incidence angle. [Ref. 6] 














Some indication of the task before the algorithm can be seen be-examining the 
frequency spectrum of that portion of-the thin-wire data-records-used by the Cadzow- 
Solomon data matrix. Figure 19 shows these spectra for. all-three aspects-available. 
Notice the varying levels of excitation of the modes. Notice-also that the spectrum 
is rather complicated considering one expects frequencies based on just integer 
multiples of the half wavelength matched to the wire length: 1.5, 3.0, 4.5,-6.0 GHz, 
etc. This is due, in part, to the fact that some ten data points at the end of early- 
‘time are included to gain the most excited resonances. This early-time data is-used 
‘by the algorithm to model the response of the feed-forward part of the target transfer 


function to the excitation waveform. 


3. Aircraft Extracted Poles 

The poles extracted from the model aircraft scattering data using the 
Cadzow-Solomon algorithm were much more highly sensitive to the input parameters 
and were less likely to group as well as the thin-wire data. The pole locations, used 
later in filtering, represent an average of the extracted location for various 
combinations of algorithm parameters and incidence angles. Individual extractions 
resulting in an outlying location from a cluster or outside the unit circle were ignored. 
The highly complex nature of the aircraft scatterer made difficult work for the 
algorithm. Figure 20 shows the frequency response of the portion of the data record 
processed for aircraft target two for the nose-on, tail-on, and thirty degree aspect. 
The spectrum is rich with closely spaced, strong frequency components at low 


frequencies which the Cadzow-Solomon algorithm must sort out. Figure 21 shows 
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Figure 18. Extracted poles, measured thin-wire data [Ref. 6] 
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this same spectral information for target three. Notice that the low frequency cuntent 
of both aircraft are similarly populated and complicated. This made the 
extraction/filtering discrimination problem very difficult, as all aircraft models tended 
to have a similar frequency spectrum, in terms of many significant components at 
lower frequencies. These lowest order poles would ideally be the ones used: for 


target discrimination later in the filtering process because they have the most energy. 


[Ref. 6] 


C. IMPLEMENTATION CONSIDERATIONS 
Actual implementation of a transient pole extraction facility must address the 
questions of frequency-subject matching and subject accuracy, signal-to-noise ratio 


required, and post-extraction pole position processing. 


1. Scattering Subject and Frequency 

The scattering objects-used for pole extraction must faithfully represent the 
actual targets of interest that will be encountered by the radar/annihilation filter 
system. This may only require more detailed scaled models than used above or may 
require full-scale targets. The two advantages of using scale models are the linear 
increase in frequency used (and linear decrease in size of antennas and components 
required) and the "reduced" availability of full scale targets — especially of aircraft of 
a potential adversary! However, scale models used must accurately represent the 
current distribution properties of the full scale targets. Validation of the modeling 


may be conducted on an aircraft available for both full-scale and model testing. 
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Figure 21. Spectra of target No. 3 data processed for pole extraction 
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Referring to the dimensions of the aircraft targets used in this thesis, Table 
4, Chapter II, indicates the low order resonant frequencies of full scale-targets lie in 
the 7 to 70 MHz region. This frequency range-would be much more difficult-to work 
in, not only due to the size of the antenna components, but also due to the 
requirement to create a broadband excitation pulse or swept frequency CW system 
in order to excite a number of the target’s resonances. The ability to create a high 
power, short pulse to-initiate the resonances is difficult at these low frequencies. The 
use of steady-state, stepped-frequency excitation may be.a way to avoid this problem 
[Ref. 25]. Here, the target is excited by relatively long pulses or CW with 
narrow bandwidth subpulses, each centered on various frequencies in the resonance 
region. The scattered response is received coherently for each frequency and 
extruded from the fixed clutter, transmitter, and receiver response by deconvolution, 
in a similar manner to that done in the TESL (Chapter II, Section B). The target’s 
broadband transient response is synthesized from the steady-state response at the 
various transmit frequencies by smoothing-in the frequency domain and subsequently 
transforming back to the time domain. This also allows tailoring the bandwidth of 
the excitation to those frequencies of most interest, vice transmitting an indiscriminate 


broadband pulse [Ref. 26]. 


2.  Signal-to-Noise Ratio 
As stated above, the pole extraction process requires approximately 20 dB 
signal-to-noise ratio. The extraction process inherently requires higher SNR than 


more conventional RF imaging from ramp responses because of reduced integration 
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opportunity; the necessity. to process mostly Jate-time, undriven responses; and the 
sensitivity of the system poles to perturbations [Ref. 27]. This requirement 
may dictate the choice between full or scale model targets. The frequency range of 
resonances of full scale targets of interest, generally 2 MHz to 1 GHz [Ref. 25], 
contains considerable man-made and atmospheric noise [Ref. 28:pp. 34.2- 
34.13], driving received noise levels well above the predominantly thermal noise 


region encountered for the resonant frequencies of scale model targets. 


3. Post-Extraction Position. Processing 

As demonstrated with the model aircraft data, the extracted pole positions 
fora given set of data vary with the processing parameters chosen. Additionally, 
although the transfer function of the target is aspect independent, the relative 
excitation of the target’s modes affects the extraction algorithm’s output. The poles 
extracted by Murphy [Ref. 6], as discussed-in the previous section, were gathered 
through a process of empirically optimizing a set of extraction parameters, then 
dithering them and averaging. The extracted pole positions-are defined by some 
probabilistic distribution, therefore the "best" set of poles for-a given target, relative 
to its actual, unknown poles, is not defined. Artificial intelligence or an expert system 
would be a logical choice to process and refine the position of extracted poles for a 
real world implementation. The expert system stage strapped to the output of the 
pole extraction algorithm could perform the sifting of extracted poles for outliers 
(performed by human eye for the poles above [Ref. 6]) and subsequent averaging. 


Possibly, a feedback path could be initiated to allow the smart processor to guide the 
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selection of extraction algorithm parameters to optimize the extracted positions. 
Additionally, this path could include the evaluation of the performance of the 


resultant annihilation filters to refine the poles. 
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IV. ANNIHILATION FILTERING 





After extracting the necessary poles for the targets of interest, using the 
algorithm of the previous chapter, these poles were incorporated into the design of 
annihilation filters. The two alternatives for this process explored here are the K- 


Pulse resonance annihilation filter and the inverse ARMA filter. 


A. RESONANCE ANNIHILATION FILTERING 
The scattered response of a target, as discussed in Chapter I, consists of the 
early and late-time portions as well as any noise collected. Morgan describes the 


transition point, T,, with unit step functions [Ref. 7], 


yt) = y(O[u()-ut-T,)]+y,Ou(t-T,)+NO. (14) 
The early-time component is composed of excitation-dependent physical optics 
scattering and natural mode startups, while the late-time portion consists of the 
natural mode resonances. This late-time portion can be written as a sum of damped 


sinusoids, 


y,@) = SOA, exp(a,t) cos(w,t+9,) (15) 
p=l 


where the amplitudes and phases (complex residues) of each mode, A,exp(/g,), are 
aspect-dependent. The additive noise and clutter component, N(f), is composed of 


atmospheric and thermal noise constituents, as well as other possible signal pollu‘ants. 
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The resonance annihilation filter (RAF) concept employs a bank of filters, 
Figure 22, with the output of each being the result of convolution of the filter impulse 


response, h,,, and the scattered response, y(t), from Equation (14):above, 


w(t) = h,(O*y(). (16) 


The decision criterion for deciding the type of target present is to match the filter 
output with the lowest average late-time energy, 

T; 

i w>()dt (17) 

2 Rae 
" Ts-T, 

with the resonances embodied in that filter. Here 7, is the start of late-time, picked 
late enough to ensure no early-time interaction remains, and 7; is the end of the 
integration time, which would ideally approach infinity for maximum discrimination, 
but must be time-limited to prevent noise from obliterating the resonant energy. This 
late-time start point is extended beyond the target late-time discussed in Chapter I, 
Section A, by the delay of the annihilation filter, 7, Here the input signal, y(t), is 
normalized by its late-time RMS energy and the filter response, H,,,, is normalized by 
its RMS value over its record length [Ref. 10]. Morgan discusses the requisite 
features of the filters, A,,, in [Ref.7] and points out that using a finite impulse 
response (FIR) filter precludes the possibility of early-time energy being carried over 


into late-time and polluting the target discrimination process. 
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Figure 22. RAF identification process [Ref. 7] 











B. THE K-PULSE FILTER 


1. Methodology 
In order to realize the annihilation filter processing with digital computers, 
the discrete time filter impulse response, h,,, must be solved for. The following 
derivation is due to Morgan [Ref. 29]. First the late-time response in 
Equation (15) is written in discrete time format, utilizing the complex z-transform 


operator to represent the sinusoidal poles, 


yn) = YO az”, (18) 


where z,=exp(s,A/), z_,=2z; =exp(s; Af), s represents each pole in the Laplace 
domain, and a, is the combination of the magnitude and residue terms which scale 
each pole’s relative contribution to the scattering. Additional flexibility can be 


incorporated by allowing k integer delays between filter weights, after which the 


convolution process of Equation (16) is written as a difference equation 


m=M 


w(n) = > h,,y,(n-kn). (19) 


m=0 


After substituting the late-time response of Equation (18) and rearranging the order 
of summation, the equation below is set to zero, indicating the annihilation action on 


the late-time energy, 
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pe-N 
p#0 


p=N ‘m= Me : 
wn) = >> zp (5: bs”) = 0. (20) 
m=0 7 


The coefficient, 1, can be set to normalize the noise power passed through the filter, 


yielding the result 


m=M 


> haz = hy. (21) 


m=} 


This result can be adjusted to represent a real linear system while reducing the 


dynamic range of the-coefficients, which results in two sets of N equations 


Men 
h,, fexp(o, 02 cost Mon}, t = -hy cost Mo, (22a) 


May 
h,, fexp(0, 9)" sin Mma = -h, sin Mac (22b) 


where p = 1 to N. To form a square matrix of equations which results in exact 
cancellation, M must equal 2N. This derivation essentially places filter zeros upon 
the poles of the input signal to be annihilated. 

Simon [Ref. 10] discusses a different approach to achieve an equivalent 
result. The coefficients, a, and £,, are solved for each mode by substituting the 


discrete time version of Equation (14) into a three-point difference equation 
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ay [(n-NAt]+y,[nAd]+ By [+I Ag = 0. (23) 
This result yields an individual three-point FIR filter for each mode. The z-transform 
derivation outlined above is equivaient to combining these in cascade and overlapping 
the leading and trailing points for each filter to form a single filter to handle all the 
target’s modes. 

The program to create the K-Pulse filter impulse response based on pole 
pair input is in Appendix D. Figure 23 contains an example K-Pulse filter impulse 
response created for the medium Q synthetic pole set; the other K-Pulse filters look 
similar. 

2. Smoothing 

The K-Pulse filter algorithm is difference equation based, which results in 
noise amplification through the filter. Notice that in Figure 24, which displays the 
spectrum of the filter of Figure 23, that the filter rises with frequency past the region 
where it performs the annihilation. In fact, the dynamic range of this.plot precludes 
examining the structure of the filter’s spectral magnitude at low frequencies, The 
frequency rise at the high frequency end of the spectrum, after the last zero, goes as 
f°", where N is the number of pole pair zeros incorporated in the filter [Ref. 29]. If 
the sampling bandwidth cannot be made small enough that the last zero is near the 
f/2 folding frequency, smoothing of the filter’s impulse response must be 


accomplished before it can be used on data with noise. Simon convolved the K-Pulse 


filter impulse response with a double Gaussian smoothing function [Ref. 10]. 
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Figure 24. Spectrum of medium Q K-Pulse filter 








Referring back to Figure 4, Chapter II, the spectrum of the double Gaussian is band- 
limited with the high and low frequency roll-offs determined by the width:of the rapid 
and narrow Gaussian pulses used in its formation. The-high frequency tail drops off 
as a Gaussian, that is exp(-af”), where a is a coefficient determined by the relative 
widths of the pulses. This decaying exponential is used to overcome the polynomial 
growth of the K-Pulse filter. The programs used to create the smoothing functions 
and the smoothed K-Pulse filter are in-Appendices E and F. The bandwidth of the 
smoothing functions used was determined by qualitatively examining the spectrum of 
the smoothed pulse. The desired filter frequency response would indicate the highest 
frequency zero with a definitive local minimum and then the spectrum would rapidly 
fall off, driven by the double Gaussian. The smoothed version of the ARMA medium 
Q synthetic data K-Pulse of Figure 23 is in Figure 25; its spectrum is Figure 26. 
Smoothed K-Pulse filters were used to test the annihilation filter concept on all the 


various signals discussed in Chapter II. 
3. Performance 


a. ARMA Synthetic Scattering 
To evaluate the effectiveness of the K-Pulse filter, scattering data was 
created for the poles sets indicated in Chapter II, as well as for poles both five 
percent above and below the baseline sets, at infinite (no noise), 20, 15, 10 and 7 dB 


signal-to-noise ratio. These data records were convolved with the K-Pulse filters 


created for the "correct" pole pairs using the convolution routine in Appendix G. As 
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Figure 25. Smoothed medium Q K-Pulse filter 
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Figure 26. Spectrum of smoothed medium Q K-Pulse filter 











a representative case, Figure 27 shows the output of the medium’Q smoothed K- 
Pulse filter operating on the-20 dB correct:pole scattering data. Figure 28:shows the: 
same filter processing the 20 dB five percent above data and Figure 29 shows the 
same filter processing the correct data at 7 dB SNR. Table 6 shows the resulting 
ratios of average late-time energies between the correct and incorrect target/filter 
matching. 


TABLE 6. RATIOS OF LATE-TIME ENERGIES: ARMA SYNTHETIC DATA 





Test Signal Signa! to SNR 
Filter 
_ Match ; * . 
re 20 dB 15 4B 10 4B 7 4B 
High Ql 5% above —-28243 11.03 4.25 2.07 156 
5% below 74442 26.80 8.53 3.59 2.29 
Medium Q? 5 % above. —-:20639 4.05 1.91 1.26 1.03° 
5% below 67721 11.79 44 2.07 1.53 
Low Q? S%above 3568 0.98 0.98 0.99 0.99 


5 % below 8814 1.03 0.99 0.99 0.99 





Note 1: Integration time = 5.14 ns 
Note 2: Integration time = 1.64 ns 


Late-time commences at T,=2.86 nanoseconds for these trials. This 
includes the 0.3 nanosecond incident pulse width after a 0.2 nanosecond delay, 
approximately 0.78 nanoseconds resulting from 20 stages in the ARMA signal 
generator which represents the 2D/c propagation delay over the scatterer, and an 
additional 1.57 nanoseconds of "shadow-time" created by the delay of the smoothed 
K-Pulse filter. The length of the integration times beyond the start of late-time was 


determined by observing the result of the mismatched filter. The summation period 
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was terminated at a point just past-where-the-significant signal structure remained in 


the output of an incorrectly matched filter to avoid noise "swamping out" 


discrimination. This was determined-to-be eight nanoseconds into the time record: 


for the high Q data and 4.5 nanoseconds for the medium and low Q case; this 
resulted in integration periods of 5.14 and 1.64 nanoseconds, respectively. Ratios of 
3 dB are still achievable down to 10:dB SNR for the medium Q five percent below 
target. This indicates the validity of:the K-Pulse filtering discrimination for a simple 
target case in the presence of noise. a 

b. TDIE Thin-Wire Scattering 

The synthetic thin-wire data was evaluated by creating a K-Pulse filter 
matched to the extracted poles, as well as-filters derived from five percent variations 
in the extracted poles. Because of the lower inherent sampling frequency of the 
TDIE program which created these signals, 480 samples in 20 nanoseconds, and the 
fact that poles were extracted at sufficiently high enough frequencies, these kill pulses 
did not require smoothing. The correct filters were convolved with both correct and 
incorrect data at various signal-to-noise ratios. The results are summarized in Table 
7. The late-time integration ran between 1.0 and 3.0 nanoseconds. 

The results of the thin-wire data point to the aspect independence of 
the annihilation filter concept. Notice that the ratios are greatest for low aspect 
angles. This is due to the fact that a single late-time start value was used, when in 
fact, for the broadside (90 degree case) there is essentially no 2D/c contribution to 


late-time. Therefore, the latest late-time value selected, which is based on the longest 
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Figure 27. Filter output for correct, medium Q synthetic data, 20 dB SNR 
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Figure 29. Filter output for five percent above, medium Q data, 7 dB SNR 














TABLE 7. RATIOS OF LATE-TIME ENERGIES: TDIE. THIN-WIRE DATA 


Target Aspect Filler-to SNR 
(Degrees) Target 
Match 
Ps 20- 1 0 4 : 
5 5% above —-214.58 11.72 1.59 216 1358 


5% below 462.50 —- 24.46 31 3.68 2.42 2 






























5% above 
5% below 


50.92 





12.52 
21.09 


5.68 
8.68 





96.37 














45 5% above 91.80 14.92 6.44 292 2.00 
5% below 132.90 19.94 8.04 3.30 — 2.13 
70 5% above 531 456 3.60 239° 1.82 
5% below 5.33 4.72 3.82 2.64 2,06 
90 5% above 56.59 336 182 4.27 1.13 
5% below 49.94 2.48 1.42 1.12 1.07 


target presentation aspect, starts much later than the significant late-time portion for 
the broadside case. Additionally, for the symmetric broadside aspect, only the odd 
modes are excited, further reducing the discrimination. 

Figure 30 contains the correct and five percent above filter outputs 
for the 30 degree aspect case at 15 dB SNR. Notice the residual energy. past one 
nanosecond for the incorrectly matched filter. Figure 31 shows the spectra for the 


late-time of these two output data records. | 


c- Measured Thin-Wire Scattering 
The same five percent variation of all-aspect, averaged poles was | 
incorporated into filters created for the measured thin-wire data. The input data 


records for all three aspects (30, 45 and 90 degrees) were passed to all three filters. . 
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Figure 30. TDIE thin-wire filter output, 30° 
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Figure 31. Spectrum of 30° TDIE thin-wire filter output 
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The resulting ratios of Jate-time energies are in Table 8. 

The best case for this data was 30 degrees aspect,.comparing the five 
percent below to the correct filter. The comparisons of the annihilation filter-outputs 
in the time and frequency domains are in Figures 32 and 33, respectively. Notice that 
even the correct filter output still contains some frequency information. This is due 
to the biasing of the extracted poles from their true location in the Cadzow-Solomon 
extraction process. These relatively small inaccuracies, however, do not preclude 
proper identification of members of this simple class of targets. 


TABLE 8. RATIOS OF LATE-TIME ENERGIES: MEASURED THIN-WIRE 








Target Aspect (Degrees) Filter to Target Match Ratio 
30 5% above 10.22 
5% below 12.68 
- 45 5% above 481 
5% below 636 
9” 5% above 271 
5% below 2.28 


zZ 


d. Aircraft Scattering 
The extracted poles for aircraft target No. 2 were incorporated into 
an annihilation filter, which was then smoothed. Scattering signals for various aspect 
angles of targets No. 2 and 3 were passed through the filter. Late-time energy 
summation was performed for two nanoseconds starting at the appropriate time, 


based on the late-time inherent in the scattering signal plus a processing "shadow- 





: time” resulting from the K-Pulse filter. The smocthed K-Pulse filter for target No. 
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Figure 32. Measured thin-wire filter output, 30° 
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Figure 33. Spectrum of 30° measured thin-wire filter output 
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2 is shown in Figure 34, its spectrum, in Figure 35. The results-of this test, in Table 
9, show that the filter could not correctly discriminate target No. 2 from target No. 
3. This problem again relates to the difficulty in extracting all the low frequency 
poles for each target and the similar high population of significant frequency 
components in this region for both targets, as-discussed in Chapter III. Figure 36 
shows the output waveforms for targets No. 2 and No. 3, each with nose-on incidence 
over the late-time energy integration period. Figure 37 shows the spectra of these 
output data records. Notice that a pole at approximately 2.5 GHz was not 
annihilated in the correct target signal. This was due to the fact that this pole was 
extracted for both targets, and, as such, was not included in the target No. 2 filter. 
At this specific aspect, however, it is excited in targec No. 2 and not target No. 3, 
which leads to incorrect identification. Figures 38 and 39 contain the time and 
frequency domain plots for the best case indicted in Table 9, the planform view with 
the wings parallel to the electric field. Here, while the 2.5 GHz resonance seems 
present in both aircraft, it does not preclude correct identification. These late-time 
only FFT plots indicate that the effective signal-to-noise ratio in late-time is much less 
than that if the entire signal, including the strong early time return, is included. 
Generally, for all the targets examined in this thesis, the ratio of early to late-time 
peak power is approximately 12 to 16 dB, therefore the late-time signal-to-noise ratio 
is less than the early-time SNR by that much. 

The observations above relate to another result encountered here, 


namely, a lack of frequency resolution in the K-Pulse filter algorithm was detected. 
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Figure 34. Smoothed K-Pulse fi 
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Figure 35. Spectrum of smoothed K-Pulse filter, target No. 2 
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Figure 36. Late-time output of targets No. 2 and No. 3 from target No. 2 
filter, nose-on 
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Figure 37. Spectra of output of targets No. 2 and No. 3 from target No. 2 
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Figure 38. Late-time output of targets No. 2 and No. 3 from target No. 2 
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Figure 39. Spectra of output of targets No. 2 and No. 3 from target No. 2 
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TABLE 9. RATIOS OF LATE-TIME ENERGIES: AIRCRAFT DATA 


Target Aspect (Degrees) Ratio (Tgt 3. LTE/Tgt-2 LTE) 





Nose On 0.20 
30 0.38 
90° 021 
Tail On 3.33 
Planform, Wings Parallel 338 
Planform, Fuselage Parallel 1.32 


These targets (targets No. 2 and No. 3) present very closely spaced poles upon 
extraction, as discussed in Chapter III, Section B. The target No. 2 filter used for the 
filtering results in Table 8 contains a zero at 1.3 and 1.5 GHz. Examination of a 
"zoomed-in" portion of the FFT of the unsmoothed version of this filter, Figure-40, 
shows that the zero corresponding to 1.5 GHz is not clearly distinguishable from that 
at 1.3:GHz. This means that any target with a pole in this neighborhood will have 
its energy annihilated, thus reducing the probability of correct discrimination. For 
this reason, the pole at approximately 2.5 GHz, separated between targets two and 
three by less than the poles discussed-here, precluded incorporation of this zero into 


the target two filter. 


e. Summary 
These results are not encouraging, however, further work in the 
refinement of the extraction process and improvements in the annihilation filters may 
enable the use of this concept on real scattering data, either from models or full-scale 


aircraft, down to near the signal-to-noise ratios realized by the synthetic data. The 


79 














x10 


Magnitude 





0 0.5 1 1.5 2 2.5 
Frequency (GHz) 


Figure 40. Close-up of FFT of target No. 2 K-Pulse filter 











aircraft targets’ frequency responses reveal them to ‘be extremely complicated. 


scatterers relative to the other, simpler models used in this thesis. ‘Concentration on 
intermediate levels of target complexity will allow these algorithms to mature or 
other, more applicable, ones to surface. One possible avenue of investigation is 
reducing the bandwidth of the excitation in order to limit the number of similar poles 
extracted and filtered for each different target. Investigation of the frequency content 
of scatterers in smaller bandwidths should be conducted to discern if there is an 
optimum frequency range where various target frequency responses offer the 


maximum potential for discrimination. 


4. lnplementatioa Considerations 
The relevant concerns for application of the resonance annihilation filter 
radar target classification method in an actual system include frequency, bandwidth, 
SNR required and its affect on transmit power, antenna design, clutter, 


electromagnetic interference, and post-filter processing and display. 


a. Frequency 
The transmit frequency for the annihilation filter radar must coincide 
with the resonances of full-scale targets - approximately 7 to 70 MHz for fighter-sized 
aircraft. The same considerations raised earlier for the scattering facility apply here 
(Chapter III, Section C), except that the option to scale the subjects to microwave 


frequencies no longer exists. The difficulties encountered at this frequency range 


include increased noise levels, the requirement for large antennas to create high gain, 


























a crowded spectrum of communications users, and difficult propagation: properties 
as the bandwidth spans the middle HF region, which includes surface and sky-wave- 
modes of propagation, to the lower VHF region, which. | is line-of-sight. 
[Ref. 30, p. 529}. . 
b. Bandwidth 
In order to effect discrimination of similarity sized targets, presumably 
with only slightly different pole locations within a-band, the system bandwidth must 


be sufficient to-excite some number of poles, possibly two to four. As stated above, 


the frequencies of interest cover a 60 MHz bandwidth. Auton indicates an ever: 


wider bandwidth of 20 MHz to 1 GHz and suggests a stepped frequency system which 
is the most economical way of creating these wide bandwidths in this frequency range 
[Ref. 25]. However, as suggested above, perhaps a subset of this-wide band would 
be an area in the frequency domain where maximum discrimination is possible among 
similarly sized targets. If so, this would ease the bandwidth requirements discussed 
here. 

Pole extraction by Murphy [Ref. 6] tends to indicate that the poles of 
differert scatterers, while different from one another, tend to reside in the same 
neighborhood for éach successive mode. While this tends to make target 
discrimination more difficult, it could allow skipping all but the most relevant 
frequencies in a stepped frequency system such as that proposed by Auton [Ref. 25]. 


This skipping would save processing time and reduce interference. The transformed 

















transjent time response could-be "filled-in" with-the missing frequency through-fitting: 


algorithms during the deconvolution process [Ref. 25]. 


Auton’s resonance based. target identification system specifically 


proposes a bank of parailel-CW transmit/receive systems each operating at'successive 


frequencies in the band of interest, each with a narrow receiver bandwidth of 
approximately 0:1 Hz. The system relies on: a companion microwave radar to 


generate range-rate information to center the narrow receiver on the target doppler 


and avoid clutter. The return signal is processed through deconvolution, as described: 


earlier (Chapter III, Section C). Auton then proposes a rational fitting algorithm to 
determine the pole locations for comparison to known locations. ‘Here, the 
substitution of a Resonance Annihilation Filter bank-may be made to speed the final 


step of the classification process. [Ref. 25] 


c. Signal-to-Noise Ratio and Range 
As a minimum, a late-time signal-to-noise ratio of 15 dB will be 
assumed for target identification via the RAF concept. Auton [Ref. 25] cites 26 dB 
at the end of an approximately ten second integration time. This has implications for 
the range equation applied to this radar. Assuming that the resonance based 
identification radar will be used with an accompanying microwave surveillance radar, 
it can be confined to operate as a surface-wave radar in the HF portion of the band 


and a line-of-sight VHF radar in that band. The range equation for this case 


(ignoring ionospheric losses) is [Ref. 30, pp. 530-531] 




















avitur 


(4n)°kTF (SIN)L, 


pi = PoP G MoT. (24) 
where R is. range, P;, is average power, G, and G,-are the t ~nsmit and receive 
antenna gains, A is the wavelength, o-is target cross section, T,-is coherent processing 
time, k is Boltzmann’s constant, F,-is received noise level above kTB due to other 
noise sources (discussed below), (S/N) is. the required signal-to-noise ratio, and _L, 
accounts for system loss. 

This frequency range, 20 to 100 MHz, is dominated:in the low-end by 
atmospheric noise and at ‘higher frequencies, by galactic noise. Assuming a rather 
large system, such as this radar identification system, would be:located ina rural area 
to avoid man-made noise, the values:of F,, in decibels above ktB, range-from 32 dB 
at 10 MHz down to 7 dB-at 90 MHz [Ref.28:p. 34.9]. Substituting these values into 
the range equation above yields Figure 41. ‘Man-made suburban or urban noise 
levels, if included-in the above calculation, would-cause the power levels required-to 
be some 20 to 30 dB higher for frequencies above 25 MHz. The other values 
involved in the equation are included in the figure and each is-discussed below. 
Nominally, such a system would require approximately 250 kW average power in 
order to sufficiently the excite the less damped poles in the low frequency end.of the 
resonance spectrum and to allow for increased: man-made noise levels-above:galactic 
noise in the frequency region above 25 MHz. More power may be required at the 


lower frequencies even still-to overcome the losses. due to surface-wave attenuation 

















caused by the radar energy coupling into the ground. This phenomenon is most 
critical around: 15 MHz and-falls off with increasing frequency [Ref. 31]. 

The radar cross section value used above Tepresents-an average; in- 
fact, as this is the resonance region for the targets of interest, the return signal 
strength can vary greatly with:small changes in frequency. As the radar generates a. 
wideband signal, the-stronger resonances, used for.identification, will boost the return 
signal. The average value indicated is derived from averaging a value at each: 
frequency determined by approximating the actual target backscatter cross section 
with that of a 2:1 prolate spheroid with perpendicular electromagnetic excitation. 


incident at 30 degrees (Ref. 32:p. 361]. 


d. Antenna Design 

The necessity to operate over a ten-to-one bandwidth will demand 
much from the antenna system used for a resonance identification radar. The-gain 
value used above, although realistic-[Ref: 30:p. 532], will not remain constant over the. 
band. As such, the-antenna will represent a:compromise design-rather than a point 
design. 

The polarization of the system would ideally be -horizontal, .as this 
would allow for the most.coupling into the resonances of the fuselage and wingspan 
of the aircraft targets in level filght. However, due to the-size of the arrays involved 
and the higher surface-wave attenuation for horizontal polarization for the lower 
frequency range of the system, the.most suitable antenna will likely be a large array 


of vertical monopoles-or dipoles, 
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In order to-reduce the clutter and noise, the radar should have: as. 


narrow azimuth and elevation beamwidths as possible. However, the-large antenna: 


aperture length of 1700-meters at 10 MHz-or 190 meters at 90-MHz. One 
alternative often exploited=at these low frequencies is the use of a broad; low gain: 
transmit beamwidth generated by driving only a-few elements of the array, overlaid: 


with a number of more precise receive beamwidths. This reduces the cost and 


complexity and’still retains:some gain-and resolution cell precision [Ref. 30:p. 531]. 


‘The varying propagation effects -bétween. the frequencies-above and: 
below 40. MHz.may warrant two separate antenna.systems:tied together for wideband. 


processing and annihilation filtering. Here, the stepped frequency nature of the 


system suggested by Auton:[Ref. 25] may make this possible. The higher-frequency 


resonances of the target scatterers.are the most highly dampéd, soit may-be found 


upon testing full-scale aircraft in actual:conditions that the the:higher end of the band 


is not cost effective. 


Another consideration-in the design and placement of the antenna is : 


the interference phenomena due to multipath propagation from reflected:signals at 
low grazing angles. The patterns which develop are highly frequency and geometry 


dependent. At the maximum extent of the effect, the return-signal strength can be 


doubled or eliminated [Ref. 33:pp. 341-361]. This effect may be ameliorated 


by the wideband nature of the system, but may. also inconsistently attenuate critical. 


resonances required for target discrimination. 
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e. Clutter 

The problem-of clutter-will-be even-more-insidious:in the resonance 
identification radar than-in normal search or tracking radar in that the frequency 
content itself of the return: signal is “mportant. Auton suggests the use of 
narrowband, doppler tracking filters to avoid miost of the clutter [Reéf. 25]. 
Additionally, controlling the propagation and precise-time-gating may reduce clutter 
from its various. sources over the band: reflections-off the surface and sky-wave 
reflection from:the ionosphere at higher HF frequencies. 


Sea. clutter is, in fact lower, in the frequency region of interest than 


in the microwave frequencies, however,. the advantage is more pronounced with | 


horizontal polarization, which:may not. be an option, considering the element. size. 
In an attempt to classify low flying: targets-over water, controlling the sea clutter by 
manipulating elevation angle may. become difficult in light of the requirement to 
avoid multipath interference. [Ref. 30:p. 483] 

The scattering data collected for this work uses a deconvolution 
process to remove non-target returns. This process may have application in the real 
world systems whereby the system could alternate the timing or direction of the 
information processed to allow a comparison of the target return and.a non-target 


return via deconvolution. This processing would have to be-done rapidly and may. 


complicate the annihilation filtering process, but may be a necessary-clutter rejection 





tactic. 
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f. Electromagnetic Interference 
Examination of a frequency allocation chart:between 7 and. 70:MHz 


indicates a crowded: spectrum of communications users filling the -band: 





F {Ref. 34:p. 1.10]. Also, due to -ionospheric effects, the generation: of ‘high: 
power, wideband: transmission :by this. system will create ‘interference problems. 
Again, a:stepped. frequency system would allow skipping protected or overcrowded 
portions of the band with subsequent deconvolution and rational function fitting filling. 
in. Additionally,-during operations-of this system,-friendly communications may be 
coordinated to avoid interference. 
As a victim of EMI, the system must deal.with-the-crowded spectrum. 
and man-made noise, as discussed earlier. This drivés up. the required: transinit 
, power for.a given SNR, and will dictate deployment location in remote, yet, ‘tactically. 
significant areas. As.a victim. of electronic countermeasures, the system’s wide 
bandwidth may offer some-inherent jam resistance. However, the criticality of the 
various frequency components used to discriminate one-target from another andthe 
requirement for a greater SNR than required for mere detection systems, may make 
the system vulnerable. Operationally, the loss of this system to jamming is, itself, a 


hostile target identifier. 


g.  Post-filter Processing and Display 





The annihilation filters operate on the decision criterion of ratios of 
late-time energy. In the presence of noise, these ratios are probablistic in nature. 


Norton conducted a rigorous analysis of the probability density function of this 
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‘process in [Ref. 4]. In order to-operate as-a-real-time {acncek or strategic-decision 
aid the target resonance-radar thust quickly-indicate its-classification-information:to 
the operator, yet it must -account for ‘the: possibilities of incorrect ‘identification-or 
inconclusive :results. This will require high speed processing of the resonance 
-annihilation-filter outputs followed by the appropriate display. This display may take 
‘the form of a-ranking ofthe most probable-identities of the target being tracked and 
the confidence level-of the identifications listed. Present-work in-artificial intelligence 
‘systems applied to-imaging systems such as-synthetic aperture radar offer a-lead in 
‘the area of real-time man-machine interface:for rapid target identification. 
h. Summary 
The implementation ofa natural resonance target identification system 
represents no greater technological risk than other radar systems that have been 
‘developed in this frequency range. The most difficult implementation task. will be 
isolating the relevant target return for‘identification processing ‘from the entirety-of 
the radar return in-time and:space, in:order-to exclude as much noise and clutter as 


possible and to concentrate the analysis on:a single target. 


€. INVERSE ARMA FILTER 
1. Methodology 
The inverse ARMA filter represents an extension of the K-Pulse filter 


concept to include canceling early-time resonance energy as well as late-time. As 
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such, it immediately surrenders-the:aspect independence property of the annihilation: 
filter discussed -earlier. 
The transfer. function-of-a‘scatterer.is represented.as a-ratio of z transform. 


polynomials 


(25) 


The inverse ARMA filter frequency response is formed by inverting this transfer 


function 


1); [1 +Q,+¢°°* €@ -N 
H,@) = SE es ne J ab 
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‘Consider the scattered response fromthe target y(t)=h,+x(t),-due to the interaction 
of the incident pulse, x(), with the-target’s impulse response, hy. In the frequency. 
domain, Y(z)=H,(z)X(z). Now, the inverse ARMA filter impulse response, /,,, is 
derived from the inverse of the target’s:transfer function, so H,, =H". The filtering 


action is described by w(t)=y(t)+*/,,, which gives this result in the frequency domain: 


W(2) = Y@)Hy(2) = Y@)H;'@ = H@X®@H, @ = X@). Gi) 
In fact, the appropriately matched-inverse-ARMA filter will-annihilate all-early and 


Jate-time.resonant energy, leaving-only a:replica of the incident-pulse. For an inverse 
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ARMA filter which is-not:matched'to the target, i.e.H,,#H;,-no Sach cancellation will, 
occur. [Ref. 35] | 

The inverse ARMA filter achieves-another :advantage over the K-Palse 
filter in that it is, itself, an ARMA filter that inherently smoothes noisy input. As 
such, it requires no additional. smoothing as the K-Pulse filter -did. [Ref. 35] The 
inverse ARMA. program is in. Appendix H. 

2. Performance 

The effectiveness of the inverse ARMA. was tested-against the synthetic 
scatterer model discussed:in Chapter Tl. This is the only available :signal for which. 
the feedtoawird: coefficients, b,, were known. For actual signals, the Cadzow- 
Solomon algorithm would have to be modified to extract-these.coefficients. Presently 
the S(gorithu ‘aiicilaciy tie siau competmadon poRion: is optimized.to determine 
the a, coefficients; the feedback coefficients which-describe the target resonances.. 
The requirement to know:both the feed-forward and feedback coefficients forces the 
inverse ARMA identification process to be aspect.dependent. 

The results of the.an inverse ARMA filter processing correctly matched, 
as well as five percent above and below, synthetic-signals is summarized in-Table 10. 
In order to allow a comparison of these results to those of the K-Pulse filter acting 
on these same targets, the late-time integration periods-were kept the:same (5.14 ns 
for the high Q, 1.64 ns for medium-and-low Q). Importantly, though, this period was 
slid forward in time to approxirnately. 0.742 nanoseconds into-the record, well into. 


early-time. Comparison of Table 10 with Table 6 indicates the advantage the inverse 
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ARMA process realizes in discrimination using-early-time energy. Figure 42 shows. 
the inverse ARMA filter output for the 26-dB-medium-Q case and:Figure 43 shows 
the output for the five percent above filter‘output. These results can be compared 
& to Figures 27 and: 28-for the K-Pulse filter. Notice that the correct: filter output is 
simply a reproduction of the incident pulse; all the target return: has been 
deconvolved out of the data. Both the high and medium Q cases show the inverse 
ARMA with significantly greater ratios of incorrect-to-correct filter late-time energies 
than ‘the normal K-Pulse filter. In fact, the inverse ARMA. performance would 
probably be better if the integration time were.not limited to the same length as the 


K-Pulse filter for direct comparison. 








TABLE 10. RATIOS OF OUTPUT ENERGIES FOR IARMA FILTER 
Test Signal Signal to SNR- 
Filler Match - P . _ 
‘  @ 20 dB i5dB «LGB 7 4B 
HighQ! 5 % above 12880 57.10 1890 «663 3.80 
5%below 14310 6612,——si28.07 792 458 
Medium.Q? S%above 4420 «-s«SB2—s—“‘<i=ésédBISC(tttiti A (8G 
S%bclow 4630 6522 2210 8.10 476 
LowQ? 5 %above 417 608 243 135 


5 % below 422° 7.56 3.25, 181 1.45 


Note 1: Integration time = 5.14 ns 
: Note 2: Integration time = 1.64 ns 


93 




















oo a are ae | . oe et ne 


i 














Time (nanoseconds) 


Figure 42. Correct IARMA filter-output at-20-dB SNR 
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Figure 43. Five percent above JARMA filter output at 20 dB SNR 




















3. Implementation Considerations 


The most significant change in:using-an inverse ARMA-filter process-for 
target idéntification, over the use-of a. K-Pulse, filter-system, is-the:fact that it is vor 
an aspect independent process. Thus the feed: forward coefficients -of the target 
transfer function, the b,’s, would have-to be determined-fora family of target aspécts 
of-most concern foreach target. The RAF concept described earlier would be much 
more-complex, as the target aspect would have to be determined. by another radar, 
then the response would have to be routed to:a filter-for that-aspect for-each possible 
target. Additionally, the number of filters required to-cover the target aspects of 
interest (say,a 60 degree sector about the-nose).would have-to be-determined.during 
the extraction process. :Perhaps the coefficients for certain discrete aspects within this. 
region could be determined and. some function: fit to them to determine those in: 
between. In any case, it is unlikely that the feed-forward portion of the target 
transfer function would be exactly matched in the [ARMA filter. This would then: 
require a closed loop, adaptive system of. filters that-could start with the coefficients 
as extracted in. the filters nearest the target aspect detected. Then, the adaptive 
controller would dither or optimize the filter coefficients in real-time-to.draw out the 
greatest ratio of energies output from the-filters. This would require much increased 
integration time on the target and an extremely fast and, perhaps, expert adaptive 
controller. [Ref. 35]: 

The concerns expressed above for implementation of the K-Pulse filter, 


such as frequency, bandwidth, and signal-to-noise ratio, remain for the inverse ARMA 
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concept. Considering the increased complexity due to the dependence on target 
aspect, -practical implementation. of the inverse ARMA filter -process for target 


identification would significantly lag that of a FIR-based, aspect independent method 


* such asthe K-Pulse filter. 














Vv. SUMMARY AND CONCLUSIONS 


A. SUMMARY 


This thesis has examined a resonance-based, aspect independent radar target 


identification system to include both major subsystems: pole extraction and 
annihilation filtering. Chapter I presented:the theoretical analysis of electromagnetic 


scattering which: gives: rise to target resonances based solely on structure and 


composition, independent of incident field aspect or polarization. Chapter II 


discussed. the four classes of signals used to-evaluate the system’s algorithms, namely: 
ARMA generated synthetic data, TDIE thin-wire data, actual measured thin-wire 
data; and data collected from scale model:aircraft. Chapter. ll presented:the theory 


of the leading candidate pole-extraction algorithm: the Cadzow-Solomon algorithm, 


results of its use on the target signals, and some initial: implications for real world. 


application. Chapter IV applied those extraction results. in resonance annihilation 
filters to evaluate target discrimination and discussed relevant areas of concern for 
actual implementation of this portion of the system. Additionally, the inverse 
ARMA, an aspect dependent resonance-based target identification scheme was 


compared to the aspect independent method. 
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B. CONCLUSIONS 

The eventual goa!-of application. of a resonance-based target identification 
system is to gain early identification of non-cooperative aircraft targets. The 
scattering of these class of targets, as modeled by the scale.aircraft used in-this thesis, 
is extremely-complex and spectrally rich, as:-compared with:that of the more simple 
thin-wire targets used herein. Intermediate levels of target: seiplenG auth as few 
element thin-wire structures should be investigated via these algorithms in-order to 
refine and-mature their application. 

The Cadzow-Solomon-pole extraction:algorithm does-adequately extract poles 
of simple, thin wire targets, both simulated-and actual. It requires a signal-to-noise 
ratio above 20 dB and its result-is sensitive-to its various input parameters, especially 
the starting point of its input data relative to the end of early-time. The over- 
estimation of the system-order and resultant ill-conditioning:of the data matrix, while 
necessary to allow for noise,:does introduce-some bias in the-extracted pole locations. 
The complexity of the frequency spectrum of the scattering from more complex 
targets, such as the model aircraft, serves as a limiting factor on the algorithm’s ability 
to extract poles. 

Implementation of a pole extraction scattering facility is feasible with the major 
caveat being accessibility to relevant targets for scattering analysis. The choice 


presented is between actual full-scale targets, which will be difficult to obtain and 


work with, considering their:size; and scale model targets, whose resonance-frequency 











thoroughly evaluated‘before implementation.. 

The K-Pulse resonance annihilation filter does effectively ‘discriminate-simple 
targets with sets of poles separated by five percent in.the-presence of:-noise down-to 
10 dB SNR. However, more complicated targets present an increasingly difficult task 
because their pole locations-are many and are closely and irregularly spaced: The 
K-Pulse filter, in its simplest-form, has a-finite limit-on the frequency resolution-of 
zeros of different targets it seeks-to discriminate. This is.due to the finite width of 
a zero null in the frequency response. of the-filter. 

Implementation of -the resonance annihilation filter radar for fighter-sized 
targets will require operation in the medium HF to low VHF bands: approximately 
7 to 100 MHz. The-design challenges here are the same as for any radar system in 
this regime and with: the additional requirement to obtain as precise a spatial and 
temporal resolution as possible in order to: isolate the -target return from noise, 
clutter, and other targets. Analysis of the range equation for a system: operating:at 
300 nautical miles indicates the requirement-for 250 kW of average transmit power. 
These requirements are within the state-of-the-art for radar systems. The greatest 
technological-emphasis required will-be in the application of advanced computing:or 
artificial intelligence to control the radar system precisely and assist the operator-in 
analysis of the-returns. 

The inverse ARMA algorithm is superior to the K-Pulse algorithm, but is more 


difficult to implement in-that it is not aspect independent. This would. require a 
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modification of the pole extraction process to allow extraction of-poles and zeros of 
the target’s:transfer function. This aspect dependence-also-adds:the requirement to 
create many filters for a single target to account for aspect. This additional 
requirement would further emphasize-high speed: processing to allow tailoring of:the 


filters to account-for the feed forward, aspect and:excitation-dependent-portion ofthe 


‘target response. 


C. AREAS FOR FURTHER RESEARCH 

The identification of non-cooperative aircraft targets may be accomplished by 
other means, such as -high-range resolution or synthetic aperture radars, or other 
target characteristics. As all these technologies mature, the relevant benefits and costs 
of each must be. compared. 

The specificalgorithms discussed earlier for resonance-based identification must 
be refined by application to more intermediate levels of target complexity. 
Additionally, study of the frequency response of.these less complicated targets, as well 
as aircraft targets, may. discern the appropriate bandwidth requirerr 2nts to allow the 


most effective target discrimination possible. 


D. UTILITY 

As the technology of the attacker becomes more advanced, effective, and 
prolific, ie. stealth, longer range precision guided munitions, etc., the-technology of 
NCTR becomes more imperative. Resonance-based target identification is sucha 


technology that holds promise for development into a viable means of NCTR. 
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APPENDIX A.. ARMA-“COEFFICIENT GENERATOR 





A. PROGRAM DESCRIPTION 
This program generates-the a, and b, coefficients of a scatterer transfer function 


given by 


(A.J) 


after reading a-data file-containing the s-plane poles for each of the N pole :pairs. 
The program also reads in the sampling inform ation to allow conversion of the input 
s-plane values to the z-plane. Additionally, the program: allows-input of a_complex. 
multiple of-each z-plane-pole pair to. allow sire ulation of the-relative amplitude and. 
phase of the transfer function: poles. 

This program was written by Capt. T. J. Murphy .and Capt. P.-C. Reddy and 


uses the subroutine POLY written by. Prof. M. A. Morgan. 


B.. PROGRAM LISTING 
PROGRAM-HCOEF 


REAL*8 T, A(0:100), B(0:2),NUM1(0:100), ATMP/0:100),NUM2(0:100)- 
REAL*8 SPOLR(30); SPOLI(30);P1, MAG(30);P!4ASE(30),C(30),D(30): 
INTEGER: NPTS, NR,-N 

COMPLEX*16 SPOL(30), ZP©t:(30). 

“CHARACTER*16- FNAME, PFNAME 


PI=3.14159265 


WRITE(*,*):’Enter filenam2 far s-plane poles ' 
READ(*,120) PFNAME 


102 











10 


18 
15 


CPEN(4,FILE=PFNAME): 
READ(1,130) NR 

WARITE(*,*)"Enter time interval’ 
READ(*,*) T 

WRITE(*,*) ’Enter number of_points’ 
READ(*,*) NPTS 

DO 10 l= 1;NR 

READ(1,*) SPOLR(I): 


_ READ(4,*):SPOLI() 


SPOL(I)= DCMPLX(SPOLR()), SPOLK(I)} 
ZPOL(l)=CDEXP(DCMPLX(T/FLOAT(NPTS .))*SPOL(|)) 
write(*,*) 'z-plane pole =’,zpol(i) 

WRITE(*,*):’Enter magnitude-of pole’, 

READ(*,*) MAG(|) 

WRITE(*, *)"Enter phase:(in degrees) of ;.0%3 J 

READ (*,*)-PHASE((): 

PHASE (!) =0.0 

MAG(I)=1.0 

CONVERT MAG-AND PHASE TO:RECTANGULAR SOORD 
C(!)=MAG(i)*COS(PHASE(|)) 

D()=MAGi(I)*SIN(PHASE (I) 

CONTINUE. 

A(0)=1.0 

A(1)=-2.0*REAL(ZPOL(1}) 

A(2)=(REAL(ZPOL(1)))**2+ (AIMAG(ZPOL(1)))**2 
NUM1(0)=2.0*C(1)- 
NUM1(1)=-2.0*(REAL(ZPOL(1))*C(1) + AIMAG(ZPOL(1))*D(#)) 


NUM1(2)=0.0 
N=2 
DO 15.1=2,NR 


NUM2(0)=2.0*C (i) 
NUM2(1)=-2.0*(REAL(ZPOL(!))*C (I) +AIMAG(ZPOL(1))*D(i)) 
NUM2(2)=0.0 

B(0)=1.0 

B(1)=-2.0*REAL(ZPCL(I)) 
B(2)=(REAL(ZPOL())))**2+ (AIMAG(ZPOL(I)))**2 
DO 17 K=0,N 

ATMP({K)=A(K) 

CONTINUE ; 

CALCULATE NEW DENOMINATOR 

CALL POLY(A,B,N) _ : 
CALCULATE NEW NUMERATOR 

N=N-2 

CALL POLY(NUMi,B,N) 

N=N-2 

CALL POLY(ATMP,NUM2,N) 

DO 18 L=0,N 

NUM1(L)=NUM1(L)+ATMP()) 

CONTINUE 

CONTINUE 
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WRITE(*;*) ’Extter name of:denominator coefficient-output:file’ 
READ(*;120) FNAME 

OPEN(3;FILE=:)° NAME) 

WRITE(3; 130)-N 

DO 20 l==1,N 

WRITE(3, *) A(j 

CONTINUE 

CLOSE(3) 

WRITE(*,*) ’Enter name Of:numerator coefficient output: file’ 
REAC(*;420) FNAME 


-OPEN(3;FILE=FNAME) 


WRITE(3;130):N 
0 21:120,N 
WAITE(3,*) NUM1(') 
CONTINUE 
CLOSE(3) 
FORMAT(E16.10) 


FORMAT(A) 
-FORMAT(I5) 


END 
SUBROLITINE POLY(A,B,N) 
Muttplying- {5(0)z**2 +.B(1)*z +-B(2)} x 


{A(0)z**N: + A(1)*2**(N-1) + A(2)*2**(N-2) + .... +-A(N)} = 
C(0)z**(N+2) + C(1)*2**(N+1) + C(2)*2**N 4 aise C(N+2): 


‘Computing C(n} coefficients. and storing in-A(n)-while 


incrementing N --> N + 2 


REAL*8°A (0:1 00),B(6:2),C(0:100): 
A(0)=1.0 

Initialize.on first call to routine 
IF(N.GE:2) GO.TO 114 

N=2 

A(i)=B(1) 

A(2)=B(2) 

GO TO 44 

C(G)=B(0)*A(0) 
C(1)=B{0)*A(1)+B(1)*A(0) 
DO.22 l=2,N. 7 
C(I) =B(0) *A(I)-+A(I-1)*B(4) +A (I-2)*B (2) 


-C(N-+1)=A(N)*B(1) #A(N-1)*B(2) 


C(N+2)=AiH)*B(2). 
N=N+2- 
DO 33.1=0,N 


AW=C() 
‘CONTINUE 


RETURN. 
END 
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APPENDIX B. ARMA SIGNAL:GENERATOR 





A. PROGRAM DESCRIPTION 
This programs generates the tine response of an ARMA system via the 
equation 


‘ksh 


k=N 
yn) = Vayn-h- + Y bxn-b 
kal k=0- 


(B.1) 


where x(1)-is the input excitatior record. a, and b, are the coefficients. determined 


by the program.in Appendix A, a'd-y(7)-is.the output-data record. 


COQNQ090000 





This program was written by Prof. M. A. Morgan. 


PROGRAM LISTING 


Program ARMA2 


Computing y[n] response for N-th order ARMA filter due to_x[n] 
input, with_coefficients: 


WR a(k) for k=1,N 
FIR b(k) for k=0,L 


by M.A. Morgan 5/23/S0 


REAL*8 x(0:2047),y(0:2047),a(1:30),b(0:32) 
CHARACTER*64 TITLE, TITL 

CHARAGTER*16 FNAME 

Entering D.E. Coefficients 

WRITE(*,*) "Enter. Filename. For Recursive.a(k):’ 
WRITE(*,*) ‘Use 0 for-FIR Fitter’ 

READ(*,100) FNAME 

IF(FNAME.EQ.’0’) GO TO 15 
OPEN(2,FILE=FNAME) 

Assuming that a(i) are-reversed polarity in file. 
READ(2,110) N 
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write(*;*):’Write k,-a(k):” 
DO-442k=1,N- 
READ(2;*)-a(k) 

a(k) =-a(k): 

write(*;*): k,a(k) 


‘CONTINUE 
CLOSE(2) 


CONTINUE : 
WRITE(*;*).’Enter-Filename for:Non-Recursive b(k):’ 
READ (*,100)-FNAME 


-OPEN(2,FILE=FNAME) 
READ (2,110)-L 


write(*,*) "Write k,-b(k):’ 
DO-22°k=0,L 


‘READ (2,*):b(k) 


write(*,*)°k,b(k). 


‘CONTINUE 


CLOSE(2): 


‘WRITE (*,*)- ’Enter Filename. for-x[n] Plot File: 


READ(*,100) FNAME 
OPEN(2,FILE=FNAME) 
READ(2,100)-TITL 


READ(2,110) NX 


READ(2,120):XQ 

READ(2,120) XQ 

write(*,*) ‘Write-m, x(m):’ 

DO:24 m=0,NX-1 

READ (2,120)-x(m) 

write(®,*)-m,x(m) 

CONTINUE 

CLOSE(2} 

Setting up.Plot File for y[n] 

WRITE(*,*) ’Enter Number of Points:-(.le. 2048)’ 
READ(*,*) NY 

WRITE(*,*) ’Enter time interval (nanoseconds)' 
READ(*,*) Tmax 

Tmin=0.0 ; ae 
WRITE(*,*) ‘Enter filename for y[n] response output’ 
READ(*,100) FNAME 

TITLE= ’ARMA Filter y[n] Response’ 
OPEN(1,;FILE=FNAME) 
WRITE(1,100) TITLE 

WRITE(1,110) NY 
WRITE(1,120)-Tmin 
WRITE(1,120) Tmax 

Initializing then iterating to form y[m] 
y(0)=b(0)*x(0) 

write(1,120) y(0). 

DO 44. m=1,NY-1 

y(m)=0.0 
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120 


Lmax=MIN(m,L) 

DO 27 k=0,Lmax 
y(m)=y(m) +b(k)*x(m-k) 
Nmax=MIN(m;N) 


‘DO 33-k=1,Nmax 


y(m)=y(m) +a(k)*y(m-k) 
write(*,*) y(m). 
WRITE(1,120) y(m) 
CONTINUE 

WRITE(*;*) ‘Thats All Folks!’ 
CEOSE(1) 

FORMAT(A) 
FORMAT(I5) 
FORMAT(E12.6) 

STOP 

END 
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APPENDIX C. NOISE POLLUTION GENERATOR 


A. PROGRAM DESCRIPTION 
This program reads:in a data file-and computes the:signal-power. After.inputing 
a seed between zero and one; it uses a routine from the FORTRAN Scientific 


Libraray called STRNUM to- generate numbers under a zero mean. uniform 


distribution. This routine is -called: 1000- times for each additive noise point 


‘corresponding to -each point of the input record, in order to. drive the random 
umber to a zero-mean -Gaussian- distribution. The variance of the random -noise 
generated-is scaled to accomodate a user input signal-to-noise ratio. 

This program is based on a signal -generator originally written by Prof..M. A. 
Morgan, modified successively by Lt..L. A. Cheeks, Lt. M: S. Simon,. and Capt. T. J. 


Murphy. 


B. PROGRAM LISTING 


PROGRAM NPLUS.FOR 


PROGRAM NOTES: 

1) GAUSSIAN NOISE ADDITION TO SYN SIGNAL OPTION 

2) PROGRAM REQUIRES 8087 MATH CO:PROCESSOR TO RUN 

3) ORIGINAL SYNGEN PROGRAM WRITTEN. IN-BASIC 
BYDR.M.A.MORGAN . 

4) STWO87FOR PROGRAM WRITTEN IN MSFORTRAN:77 
BY.L. A. CHEEKS 

5) UPDATE AUG89.BY LT.M.S. SIMON 

6) UPDATE. MARQO-BY CAPT T, J::MURPHY 


QAQAIQA9Q9QAAIANNO 


***DECEARATIONS*** 
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REAL T9,T0,A,C9,P1,XOSUM,GSUM,GSUMPRM,R,CDB,X0(1024),RO(100);10(10 
+0),A0(100);P0(100),G(1024),X1 (1024),P,TZ/0.0/,DT1,DT,L2,X02SUM,X02 
+(1024),GADJ,GMEAN- 

‘ INTEGER N,Z,19,11, NCOUNT 
CHARACTER:FN*16;TITLE*16,ANS*1;HEADER*64 
DATA P¥/3.1415927/ 


***MAIN BODY OF PROGRAM*** 


***DATA INPUT ROUTINE*** 


yoaaggg 


0 WRITE(*,*) ‘ENTER NAME OF DATA FILE TO:BE CORRUPTED’ 
READ (*,280) TITLE 
OPEN(1,FILE=TITLE). 
READ (1,260): HEADER 
READ (1,240):NPTS. 
READ- (1,250) TBEGIN 
READ (1,250): TEND: 


Cc 
DO 190 1=1,NPTS 
READ. (1,250)-XO(I) 
490 CONTINUE. 
CLOSE(1) 
C 
WRITE(*,*)'DONE READING FILE’ 
. C 
XOSUM=0.0 
DO 80 I=1,NPTS 
‘ XOSUM=XOSUM + (X0(I)**2) 


CONTINUE 





XOSUM=X0SUM/NPTS 


80 
Cc 
C  XOSUM/NPTS IS AVERAGE POWER 
GC 
C  ***ROUTINE TO ADD GAUSSIAN NOISE TO SIGNAL*** 
Cc 
GSUM=0.0 
GADJ=0.0 
GMEAN=0.0 
120 WAITE(*,*)"ENTER SEED FOR RANDOM NUMBER GENERATION FOR.NOISE’ 
WRITE(*,*)"(0 < R < 1.0) 0.21289 GENERATES A.GOOD DISTRIBUTION’ 
READ(*,250) R 
DO 130-1=1,NPTS 
CALL CLMT(R,P) 
G()=P 
GADJ=GADJ+G()) 
430 CONTINUE 
C ADJUST DISTRIBUTION SO MEAN IS 0 
GMEAN=GADJ/FLOAT(NPTS) 
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156 


160 
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DO-135 |=1,NPTS 

G(\)=G(1)-GMEAN 

GSUM=GSUM++ (G(|)**2)- 

CONTINUE 
GSUM/NPTS IS-AVERAGE POWER OF GAUSSIAN-SIGNAL 
GSUM=GSUM/NPTS 

WRITE(*,*) ENTER: DESIRED SIGNAL TO’NOISE RATIO IN-DB” 
READ(*,*)-DB 

GSUMPRM=X0SUM/(10.0**(DB/10.0)) 
CDB=10.0*ALOG10(XOSUM/GSUMPRM) 
WAITE(*,*):’CHECK:DB.=',CDB 
RATIO=SQRT(GSUMPRM/GSUM) 

DO 140 1=1,NPTS 

G(I)=RATIO*G() 

CONTINUE 

DO 150 I=1,NPTS 

XO(I)=X0(1)+G(\) 

CONTINUE 


***ROUTINE TO NORMALIZE SYN GEN SIGNAL*** 


X02SUM=0.0 
DO 156 I=1,NPTS 
X02(I)=XO(I)**2 
X02SUM=X02SUM-+X02() 


‘CONTINUE 


L2=SQRT(X02SUM): 
DO 160 1=1,NPTS 
XO(I)=XO()/L2 
CONTINUE 


***ROUTINE TO STORE SYN GEN SIGNAL DATA IN DATA.FILE*** 


WRITE(*,*):"ENTER NAME OF-DATA FILE TO*CONTAIN CORRUPTED SIGNAL’ 
READ :(*,280)-FN 

WRITE(*,*)"STORING CORRUPTED SIGNAL DATA IN DATA FILE’ 
OPEN(1,FILE=FN) 

WRITE(1,260) HEADER 

WRITE(1,240) NPTS 

WRITE(1,250)- TBEGIN. 

WRITE(1,250) TEND 


DO 17:1=1,NPTS 
WRITE(1,250) X0(I) 
CONTINUE 
CLOSE(1) 


-FORMAT(IS) 
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-FORMAT(E12.6). 
‘FORMAT(A64) 
-FORMAT(A8) 
-FORMAT(A16) 


WRITE(*,*) *CHECK NEW FILE(S)’ 


‘GO TO 170 


STOP 


‘END 


***SUBROUTINES*** 
SUBROUTINE CLMT(R,P) 


THIS‘SUBROUTINE IS BASED ON THE CENTRAL.LIMIT THEOREM AND 


“GENERATES A SEQUENCE OF INDEPENDANT NUMBERS WHICH ARE 
‘NORMALLY DISTRIBUTED:ABOUT THE INTERVAL 0 TO 1 


REAL P,R 


“P=0.0 


DO 170 Il=1,1000 
CALL STRNUM(R) 
P=P+R 
CONTINUE 
P=P/1000.0 
RETURN 

END 


SUBROUTINE STRNUM(R) 


THIS SUBROUTINE GENERATES A SEQUENCE OF-NUMBERS WHICH:ARE 
RANDOMLY AND UNIFORMLY DISTRIBUTED-OVER THE UNIT INTERVAL. 
REF:-P, 280 OF FORTRAN: Scientific Subroutine Library 


NOTES: 
1) R = 0.21289 GENERATES A GOOD UNIFORM-DISTRIBUTION 


BB=1.0 
P1=R*317.0 
R=AMOD(P1,BB) 


-RETURN 


END 
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APPENDIX D. K-PULSE GENERATOR 


‘A. ‘PROGRAM DESCRIPTION 

This.program creates an K-Pulse filter impulse response after reading a file 
containing the s-plane poles to be annihilated. This program allows scaling of the 
filter impulse response as-well:as integer delays bewteen filter weights, each:to allow 
additional-degrees of freedom.in-the filter design. In this thesis, these parameters 
were set to one. 

This program was written by Prof. M. A..Morgan, with susequent modifications 


by Lt..M. S. Simon and Capt. P. C. Reddy. 


B. PROGRAM LISTING 


‘PROGRAM:KP 


PROGRAM NOTES 
1):PROGRAM-GENERATES KILL PULSES 
2):KEYBOARD INPUTS: POLES/REAL,IMAG 
3). OUTPUT TO DATA FILE: K:PULSE TIME SERIES 
4) PROGRAM: REQUIRES 8087 MATH CO:PROCESSOR: 
5). PROGRAM WRITTEN IN MSFORTRAN-BY PROF. M.A. MORGAN 
6) INTERACTIVE MODS BY LT M.S, SIMON 
7) ADDITONAL MOS BY:CAPT-P.C.:-REDDY 


K-Pulse Generation Using Z-Transform Product Non-Recursive 
Fitter Design Approach by M.A. Morgan 8/24/89, 


on0000 00000000 


-REAL-SR(20),S1(20),A(0:100),B(2),H(401),SF 
CHARACTER FNAME*16,TITLE*64,ANS*1,PENAME*16 
WRITE(,*)-'*** K-PULSE GENERATOR *** 
WRITE(*;*) . 

1 ‘WRITE(*,*)"ENTER NUMBER OF POLE. PAIRS (<=80):’ 
WRITE(*,*).'ENTER SCALE FACTOR’ 

READ(*,*):SF 
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12 
13 


ZO0O 


‘WRITE(*,*) ENTER TIME WINDOW (NANOSEC):’ 

-READ(*,*) TO: 

WRITE(*,*) ENTER NUMBER OF TIME POINTS: 

READ(*,*) NT 

‘DT=TO0/(NT-1.0) a ue 
WRITE(*;*) ENTER TIME-STEPS BETWEEN FILTER WEIGHTS (>=1):’ 
-READ(*,*) KS: 


-OPEN(2,FILE=PFNAME) 
READ(2,*) NP 

‘DO 11 N=1,NP 
‘READ(2;*) SRN) 


‘INPUT DATA-CHECK 


WRITE(*,*)'SCALE FACTOR |S ',SF 
WRITE(*,*)' TIME WINDOW.IS ',TO,’ NANOSEC: 


‘WRITE(*,*)'DATA CORRECT’ ~ 





WRITE(*;*) ENTER FILENAME OF SPLANE POLES’ 
READ(*,100) PENAME 


READ(2,*) SI(N) 
CONTINUE 
CLOSE(2) 


WRITE(*,*) °* 
WRITE(*,*) INPUT DATA CHECK’ 
WRITE(*,*)'NUMBER OF POLE PAIRS IS ’,NP 


WRITE(*;*)'NUMBER OF TIME POINTS.IS”,NT ae 
WRITE(*,*)'NUMBER OF TIME STEPS BETWEEN FILTER WEIGHTS IS ’,KS 
WRITE(*,*)’ ’ 

DO 12 N=1,NP 

WRITE(*,*)'SIGMA:’,N,” = ’,SR(N) 

WRITE(*,*)‘OMEGA ’,N,’ = ‘,SI(N) 

PAUSE 

CONTINUE 

WRITE(*,*)'IS DATA CORRECT?’ 

READ(*,200)ANS 

IF((ANS.EQ.’Y’).OR.(ANS.EQ.’y’)) THEN 





GO TO 14 

ELSEIF ((ANS.EQ.'N’).OR.(ANS.EQ.’n’)) THEN 
WRITE(*,*)’DATA ERRORS. START AGAIN’ 
GO TO 1 
ELSE ‘| 
WRITE(*,*)'TRY AGAIN! Y OR N’ 
GO TO 13 

ENDIF 





SCALING POLES 


DO 15 N=1,NP 
SR(N)=SR(N)/SF 
SI(N)=SI(N)/SF 








CONTINUE 


CALCULATIONS 


0.09 &, 


MP=0 - 
DO_22 N=1,NP 
E=EXP(KS*SR(N)*DT) 
C=COS(KS*SI(N)*DT). ‘ 
B(1)=-2.0*E*C 
B()=E*E 
CALL POLY(A,B,MP) 

22 CONTINUE 

WAITE(**) 

WRITE(*;*) "UN-NORMALIZED K-PULSE FILTER WEIGHTS: 

WRITE(*,*) 

DO:33 M=0,MP 

WRITE(#,*) M,A(M) 

CONTINUE 
PAUSE 


NORMALIZE K-PULSE TO UNIT- ENERGY 


QO0) 8. 


ENERGY=A(0)**2 

M1=KS*MP+1 

DO 44 M=1,MP 

ENERGY=ENERGY+A(M)**2 

RMS=SQAT(ENERGY) * 


DEFINING UNIT SAMPLE RESPONSE 


Q09 8B 


H(1)=A(0)/RMS 

I=1. 

DO 55 N=1,MP 

M=0- | 
50. M=M+4 | 
I=1+1 | 
IF(M.GE.KS) GOTO 55 | 
H()=0.0 | 
GO TO 50 
H(I)=A(N)/RMS 





STORING K-PULSE IN DATA FILE 


Qd9008 


TS=0.0 

TF=KS*MP*DT 

WRITE(*,*) STORING K-PULSE IN DATA FILE’ 
WRITE(*,*) : 

WRITE(*,*) ENTER OUTPUT FILE (<= 16 CHAR): 
READ(*,100) FNAME 
WRITE(*,*) ENTER HEADER FOR OUTPUT FILE (<= 64 Char):’ : 
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READ(*,100) TITLE 


OPEN(1,FILE=FNAME) 


WRITE(1,100). TITLE 


WRITE(1,110):M1 
“WRITE(1,120):TS 


WRITE(1,120) TF 


‘DO-66 N=1,M1 
‘WRITE(1 ,420):H(N) 


CONTINUE . . 
WRITE(*;*)'K-PULSE GENERATOR PROGRAM COMPLETE’ 


‘WRITE (*,*)’CHECK NEW FILE’ 


-FORMAT(A) 


-FORMAT(I5) 


FORMAT(E12.6) 


- -FORMAT(A1) 
‘STOP 


END 


‘SUBROUTINE: POLY(A;B,N) 


‘MULTIPLYING {z**2 + B(1)*z + B(2)} x 


{2**N + A(i)*2**(N-1) + A(2)*2**(N-2) +... + A(N)} = 
2**(N-+2) + C(1)*2**(N41)-+ C(2)*2**N +... + C(N+2) 


‘COMPUTING C(N).COEFFICIENTS AND STORING INTO A(N) WHILE 
‘INCREMENTING N --> N + 2 


REAL A(0:100),8(2),C(0:100) 
A(0)=1.0 

INITIALIZE ON FIRST CALL TO ROUTINE 
IF(N.GE.2) GO TO 11 

N=2 

A(i)=B(1) 

A(2)=B(2) 

GO TO 44 

C(1)=A(1)+B(1) 

DO 22 l=2,N 

C(I) =A(l) +A(I-1)*B(1)+-A(1-2)*B (2) 
C(N+1)=A(N)*B(1)+A(N-1)*B(2) 
C(N+2)=A(N)*B(2) 

N=N+2 

DO 33 I=1,N 

A()=C\()) 

CONTINUE 

RETURN 

END 

















APPENDIX E.. DOUBLE:GAUSSIAN SMOOTHING FUNCTION GENERATOR: 


A. PROGRAM:DESCRIPTION 

This program creates-a double-Gaussian-smoothing function used to-allow the 
use of the difference equation based: K-Pulse-filter in-the presence the noise. The 
program -uses input values-of the ten percent ‘height of the low and high frequency 
ends of the double-Gaussian: frequency response to determine thé:pulse widths in the 
time domain. The user defines-a threshold value to truncate the-tails of the double 
Gaussian function to limit its temporal extent. The program. uses the PLTSUB 
subroutine written by Prof. M. A. Morgan, which uses the -Grafmatics FORTRAN 
graphics -package, to create:a plot of-the result. 

This program was originally written by Prof. M. A. Morgan and was modified: 


by Capt. P. C. Reddy. 


B. PROGRAM LISTING 


PROGRAM DGSMF 


Adaptive Double-Gaussian Smoothing Function Generator. 
By M.A. Morgan 2/22/90;:Modified by P.C. Reddy 8/1/90 


2000 


REAL Y1(1024),DT,T1,T2,F1,F2,TMIN, TMAX,TMAX1, THR 
REAL*8 Y(1024),E1,E2,F(1024),A1,A2,ALFA1,ALFA2,C1,C2 
INTEGER NPTS, N 

CHARACTER*1 YN ; 

CHARACTER*64 FNAME, TITLE 

WRITE(*,*) ‘Enter No. Time Points ’ 

READ(*,*) NPTS 

TMIN=0.0 

WRITE(*,*) 
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WRITE(*;*) ‘Enter Time Window in NScc:. ’ 
READ(*,*) TMAX 

DT=TMAX/FLOAT(NPTS-1) 

CONTINUE 6 4: 
WRITE(*;*) ‘Enter Upper10%-Freq Rolloff:in GHz ° 
READ(*,*) F1. 

T1=1.47/F1 ; 

WRITE(*;*) 'Narrow:10%:Pulse:Width in NSec:= ',T1 
WRITE(*,*) ‘Enter Lower-10%.Freq Rolloff-in GHz ° 
READ(*,*) F2 

WRITE(*;*) 

T2=.31/F2 ; 

WRITE(*,*) ‘Wide 10% Pulse Width in NSec:= ',T2 
WRITE(*;*) 

WRITE(*,*) ‘Are These Pulse Widths OK ?-(Y/N) ’ 
READ(*,100). YN 

IF(YN.EQ.’N’.OR.YN:EG.'n') GO TO-44 

WRITE(*,*) 'Enter THRESHOLD ’ 

READ(*,*) THR 


ALFA1=4*DLOG(10)/T1**2 


ALFA2=4*DLOG(10)/T2**2 
A1=DSQRT(ALFA1)/(DSQRT(ALFA1)-DSQRT(ALFA2)) 
A2=Ai-1.0 

1=0 

DOJ0K=1,NPTS | 
T=(K-1.0-(NPTS/2))*DT 
C1=ALFA1*T**2 
C2=ALFA2*T**2 
IF(C1.GT.THR)THEN 
E1=0.0 

ELSE 

E1=A1*EXP(-C1) 
ENDIF 
IF(C2.GT.THR)THEN 
F(K)=0.0 

GO TO 10 

ELSE 

E2=A2*EXP(-C2) 
ENDIF 

F(K)=E1-E2 

IF (F(K).NE.0.0) THEN 
l=I+1 

Y(Q=F(K) 
Y1()}=SNGL(¥(\)) 

M=l 


ENDIF 
CONTINUE 
TMAX1=DT*N 
EGA Assumed: 
NS=1 
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Cc Impact Printer Assumed: 
NFD=1 


TITLE=’DOUBLE-GAUSSIAN SMOOTHING: FUNCTION’ 
CALL PLTSUB(TITLE,N,TMIN, TMAX1,Y1,NS,NFD) 
WRITE(*; *) "Do. You Want to Storée-These-Results-?-(Y/N) * 
READ(*,100) YN 
-IF(YN.EQ;’N’.OR.YN.EQ.’n’): GO TO 99 
WRITE(*;*) ‘Enter-File Name for Plot Data:Storage: ’ 
READ(*,100) FNAME 
OPEN(3,FILE=FNAME) 
WRITE(3,100): TITLE 
WRITE(3,110) N 
WRITE(3,120) TMIN 
WRITE(3,120) TMAX1 
-DO 88 I=15N 
-WRITE(3,120) Y(I) 
88 CONTINUE 
CLOSE(3): 
99 -CONTINUE 
WRITE(*,*) Program Completed’ 
100 FORMAT(A). 
110 FORMAT(I5) 
120 FORMAT(E12.6) 
STOP: 
END 
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APPENDIX F. SMOOTHED:K-PULSE GENERATOR 


A. PROGRAM DESCRIPTION 

This-program is an adapataion of a linear convolution program written by Prof. 
M. A. Morgan. This program reads. two data files and convolves them and 
normalizes the-output waveform energy before depositing the-result in user specified 


data file. 


B. PROGRAM LISTING 





PROGRAM CVN 


Convolution of PLOT Format Data:Files: X(n)*¥(n):to:Produce 
2(n). Provides Local Plotting and PLOT Output File. — 
By M. Morgan, 8/21/89, Modified by P. Reddy, 8/25/90. 


aOaqgng0n0 


REAL X(1024),Y(1024),Z(2048).7.1/2048), TLATE 
, CHARACTER*1 ANS, BELL 

CHARACTER*8.FNX,FNY,FNZ 

CHARACTER*64 TITX, TITY, 7iTZ 


BELL=CHAR(7) 





DATA FILE INPUT 


5EQ0 9 


CONTINUE 
WRITE(*,*) "THIS PROGRAM CONVOLVES X(n) * Y(n) = Z(n)’ 
WRITE(#,*) ‘AND PLACES 2(n) IN USER SPECIFIED OUTPUT FILE’ 
WRITE(*,*) 

WRITE(*,100). BELL 

WRITE(*,*) "ENTER NAME OF X(n) INPUT FILE” 

WRITE(*,*) 

READ(*,100) FNX 

WRITE(*,*) 

WRITE(*,*) ’X(n) INPUT FILE IS *,FNX 

WRITE(*;*) 

WRITE(*,100) BELL 

. WRITE(*,*) ‘ENTER NAME: OF Y(n) INPUT FILE” 
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WRITE (*,*) 


- READ(*,100) ENY 


WRITE(*,*) 


WRITE(*,*) 'Y(n) INPUT-FILE IS";FNY 


WRITE (*,*) 
READING INPUT FILES 


OPEN(1,FILE=FNX) 
OPEN(2,FILE=FNY) 
READ(1,100) TITX 
READ(1,120) NX 
READ(1,130) TX1 
READ(1,130) TX2 
DO-20:1=1,NX 
READ(1,130) X(I) 
CONTINUE 

WRITE (*,*) °X(t) FILE READ’ 
READ(2,100) TITY 
READ(2,120) NY 
READ(2,130) TY 
READ(2,130) TY2 

DO 30:1=1,NY 
READ(2,130) Y(\) 
CONTINUE 
WRITE(*,*) 'Y(t) FILE READ’ 
WRITE(é,*) 

CLOSE(1) 

CLOSE(2) 


CONVOLUTION ALGORITHM WITHOUT ZERO PADDING 
Written by-M.A. Morgan-5/19/89- 


WRITE(*,*) "CONVOLUTION IN PROGRESS’ 


WRITE(*,*) 

DO 50: N=1,NY 

Z(N)=0.0 

M2=MINO(N,NX) 

DO 40 M=1,M2 
Z(N)=Z(N)-+X(M)*Y(N-M-+1) 
21(N)=Z(N) 

CONTINUE 

Ni=NY+1 

N2=NX+NY-1 

DO 70:N=N1,N2 

Z(N)=0.0 

Mi=N-+1-NY 
M2=MINO(N,NX) 

DO 60-M=M1,M2 
Z(N)=Z(N)-+X(M)*Y(N-M+1) 
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CONTINUE 


TZi=TX1+TY1 
TZ2=TX2+TY2 

ENERGY=0.0 

DOSON=1N2 
ENERGY=ENERGY +2Z(N)*Z(N) 
DO 92. N=1,N2 
Z(N)=Z(N)/SQRT(ENERGY) 


ZA(N)=Z(N). 


WRITE(*,*). 
PAUSE 


*** Plotting-Z(n) *** ; 
TITZ=:Convolution Z(n)=X(n)*Y(n)’ 
CALL PLTSUB(TITZ,N2,12Z1,TZ2,TZ1,TZ2;21,1,1,0) 


WRITE(*,100) BELL 

WRITE(*,*) “ENTER NAME OF Z(n). FILE’ 
WRITE(*,*) 

READ(*,100) FNZ 

WRITE(*,*) 

WRITE(*,*) OUTPUT FILE IS ’,FNZ 
OPEN(3,FILE=FNZ) 

WRITE(*,*). 

WRITE(*,100)-BELL ; 
WRITE(*;*) ‘ENTER HEADER TITLE FOR Z(n)-FILE’ 
WRITE(*,*) 

READ(*,100) TITZ 

WRITE(*,*) - 
WRITE(*,*) ‘SAVING Z(t) TO OUTPUT DATA FILE’ 
WRITE(3,100) TITZ 

WRITE(3,120) N2 

WRITE(3,130) TZ1 

WRITE(3,130) Tz2 

DO 80 I=1,N2 

WRITE(3,130) 2(I) 

CONTINUE 

CLOSE(3) 

QUERY TO DO ANOTHER RUN 


CONTINUE 
WRITE(*,*) 
WRITE(*,100) BELL 





WRITE(*,*) 'DO YOU WISH TO CONVOLVE ADDITIONAL WAVEFORMS?’ 


WRITE(*,*) 

READ(*,100) ANS 

IF ((ANS.EQ.'Y’).OR.(ANS.EQ.’y’)) THEN 

GO TO 10 

ELSEIF ((ANS.EQ.’N’).OR.(ANS.EQ.’'n')) THEN 














99 


100- 
120: 
130- 
STOP 
END: 


GO-TO 99 

ELSE 

WRITE(*;*) "TRY AGAINE*Y" OR "NY? 
GO TO 95 

ENDIF 

CONTINUE 

WRITE(*,100):BELL 


WRITE(*,*) CONVOLUTION PROGRAM IS OVER. CHECK NEW FILE(S)’ 


FORMAT(A) 
FORMAT(I5) 
FORMAT(E12.3) 
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APPENDIX G: K-PULSE CONVOLUTION PROGRAM 


PROGRAM DESCRIPTION 


This program is an adapataion of a linear-convolution program wriiten by Prof. 


‘M. A. Morgan. This program reads two data files and convolves them and calculates 


the average energy over a specified time.of the output record. The program plots 


the result-and deposits-it in a user specified data file. This program was modified by 


‘Capt. P. C. Reddy. 
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PROGRAM LISTING 


PROGRAM CV2 


Convolution of PLOT Format Data Files X(n)*Y(n) to Produce 
Z(n). Provides Local Plotting and PLOT Output File Along 
With Late-Time Energy Calculation. 

By M: Morgan, 8/21/89. Modified by P. Reddy, 8/25/90. 


REAL X(1024),Y(1024),Z(2048),Z1 (2048), TLATE 
CHARACTER*1 ANS, BELL 

CHARACTER?*8 FNX,FNY,FNZ 
CHARACTER*64 TITX,TITY, TITZ 


BELL=CHAR(7) 
DATA FILE INPUT 


CONTINUE 

WRITE(*,*) "THIS PROGRAM CONVOLVES X(n) * Y(n) = Z(n)’ 
WRITE(*,*) 'AND PLACES Z(n) IN USER SPECIFIED OUTPUT FILE’ 
WRITE(*,*) 

WRITE(*,100) BELL 

WRITE(*,*) "ENTER NAME OF X(n) INPUT FILE’ 

WRITE(*,*) 

READ(*,100)-FNX 

WRITE(*,*) 

WRITE(*,*) 'X(n) INPUT FILE IS ',FNX 














WRITE(*,*) 

WRITE(*,100) BELL 

‘WRITE(*,*) ‘ENTER NAME-OF Y(n) INPUT FILE’ 
WRITE(*,4) 

READ(*,100) FNY 

“WRITE(*,*) 

WRITE(*,*) ’Y(n) INPUT FILE IS \FNY 
‘WRITE(*,*) 





READING INPUT FILES 


aqa0o”0 


OPEN(1,FILE=FNX) 
OPEN(2,FILE=FNY) 
‘READ(1,100) TITX 
READ(1,120) NX 
READ(1,130) TX1 
READ (1,130) TX2 
DO 20 I=1,NX 
“READ(1,130) X(I) 

20 CONTINUE 
WRITE(*,*) 'X(t) FILE READ’ 
READ(2,100) TITY 
READ(2,120) NY 
READ(2,130) TY1 
READ(2,130) TY2 
DO 30 I=4,NY 
READ(2,130) Y(I) : 

30 CONTINUE 
WRITE(*,*) ’Y(t) FILE READ’ 
WRITE(*,*) i 
CLOSE(1) 

CLOSE(2) 


CONVOLUTION ALGORITHM WITHOUT ZERO PADDING 
Written by-M.A. Morgan 5/19/89 


aagg 


WRITE(*,*) CONVOLUTION IN PROGRESS’ 
WRITE(*,*) 
DO 50 N=1,NY 
Z(N)=0.0 
M2=MINO(N,NX) 
DO 40 M=1,M2 
Z(N)=Z(N)+X(M)*Y(N-M+1) 
40  Z4(N)=Z(N) 
50 CONTINUE 
N1=NY+1 A 
N2=NX+NY-1 
DO 70 N=N1,N2 
Z(N)=0.0 
M1=N-+1-NY 














M2=MINO(N,NX) 
DO-60.M=M1,M2 
Z(N)=Z(N)+X(M)*¥(N-M+1) 
Z4(N)=Z(N) 

CONTINUE 


CALCULATING LATE-TIME ENERGY 


TZ1=TX1+TY4 

TZ2=TX2+TY2 

WRITE (*,*) . 

WRITE(*,*) ENTER TL (Nsec) FOR LATE-TIME ENERGY CALCULATION’ 
WAITE(*,*) 

READ(*,*) TLATE . 
WRITE(*,*): "ENTER TZ3 (Nsec) FOR END OF’ 
WRITE(*,*) LATE-TIME ENERGY CALCULATION’ 
‘WRITE(*,*) 

READ(*,*) TZ3 
‘NLATE=NINT(TLATE*N2/(TZ2-TZ1)) +1 
N3=NINT(TZ3*N2/(TZ2-TZ1)) +1 
IF(NLATE.LT.N2) GO TO 85 

WRITE(*,100) BELL . 
WRITE(*;*):"TLATE TOO LARGE - TRY AGAIN’ 
GO-TO 83 

CONTINUE 

IF(NLATE.GE.1) GO TO 87 

WRITE(*,1G0) BELL 

WRITE(*,*) "TLATE TOO SMALL - TRY AGAIN’ 
GO TO 83 

CONTINUE 

ENERGY=0.0 

DO 90 N=NLATE,N3 
ENERGY=ENERGY+2Z(N)*Z(N) 
ENERGY=ENERGY/(N3-NLATE+1) 
WRITE(*,*) 

WRITE(*,*) ‘LATE TIME ENERGY =’, ENERGY 
PAUSE 


*** Plotting Z(n) *** 
TITZ=’Convolution Z(n)=X(n)*Y(n)’ 
CALL PLTSUB(TITZ,N2,TZ1,1Z2,T2Z1,TZ2,2Z1,1,1,0) 


WRITE(*,100) BELL 

WRITE(*,*) ‘ENTER NAME OF Z(n) FILE’ 
WRITE(*,*) 

READ(*,100) FNZ 

WRITE (*,*) 

WRITE(*,*) ‘OUTPUT FILE IS ’,FNZ 
OPEN(3,FILE=FNZ) 

WRITE(*,*) 


~ 125 











99 


100 
120 
130 








WRITE(*,100) BELL . 
WRITE(*,*) ENTER HEADER-TITLE FOR Z(n) FILE” 


‘WRITE(*,*) 


READ(*;100) TITZ 

WRITE (*,*) 

WRITE(#,*) ‘SAVING Z(t) TO OUTPUT: DATA FILE’ 
WRITE(3,100) TITZ 

WRITE(3,120):.N2 

WRITE(3,130) TZ1 


WRITE(3,130) TZ2 


DO 80 1=1,N2 

WRITE(3,130) Z()) 

CONTINUE 

CLOSE(3) 

QUERY TO DO ANOTHER RUN 


CONTINUE 

WRITE(*,*) 

WRITE(*,100) BELL 
WRITE(*,*) 'DO'YOU WISH TO CONVOLVE ADDITIONAL WAVEFORMS?’ 
WRITE(*,*) 

READ(*,100) ANS 

IF((ANS.EQ.'Y’):OR. (ANS.EQ.’y’)) THEN 

GO TO-10 

ELSEIF ((ANS.EQ.’N’).OR.(ANS.EQ.'n’)) THEN- 

GO TO 99 

ELSE 

WRITE(*,*) ‘TRY-AGAIN! *Y! OR 'Nt' 

GO TO 95 

ENDIF 

CONTINUE 

WRITE(*,100) BELL 

WRITE(*,*) ‘CONVOLUTION PROGRAM IS OVER. CHECK NEW FILE(S)’ 
FORMAT(A) 

FORMAT(I5) 

FORMAT(E12.6) 

STOP 


END 











APPENDIX H. INVERSE ARMA PROGRAM 


‘ A. PROGRAM DESCRIPTION 





This program implements an inverse ARMA aspect dependent annihilation 
filter. As such, it is an adaptation of the ARMA signal generator in Appendix B. It 
requires input of the a, and. b, coefficients for the target transfer function it is 
intended to annihilate. It then creates the inverse impulse response and convolves 
it with an input data record. Finally, it computes the average energy over a user 
specified time period in the output. 

This program was written by Prof. M. A. Morgan and modified by Capt. P.-C. 

F Reddy. 


B. PROGRAM LISTING 


Program IARMA 
Reading N-th order ARMA filter with coefficients: 


IR a(k) for k=1,N 
FIR b(k) for K=0,L 


Processing y[n] using the inverse ARMA filter with output 
z[n]. Inverse ARMA has coefficients: 


IR ai(k)=-b(k)/b(0) for k=1,L 
FIR. bi(k)=-a(k)/b(0) for k=1,N 
bi(0)= 1.0/b(0) 





by M.A. Morgan 5/24/90 
Modified by P.C. Reddy 5/27/90 


QANDNANQNNQ0NN0NNN9NN9NN0 


y REAL*8 y(0:2047),z(0:2047) 
REAL*8 a(1:30),b(0:32),ai(1:30),bi(0:30) 
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24 





CHARACTER*64 TITLE, TITL 
CHARACTER*16 FNAME 
Entering.D.E. Coefficients 
WRITE(*;*) ‘Enter Filename For Recursive a(k):’ 
WRITE(*;*) ‘Use 0 for FIR Filter’ 
READ(*,100) FNAME 
IF(FNAME.EQ:0’) GO TO-15- 
OPEN(2,FILE=FNAME) 
Assuming that a(k) are reversed polarity in file. 
READ(2,110)-N 

write(*,*) Write k, a(k):’ 

DO 11-k=1,N 

READ(2,*) a(k) 

a(k)=-a(k) 

write(*,*) k,a(k) 

CONTINUE 

CLOSE(2) 

CONTINUE 

WRITE(*;*) ‘Enter Filename for Non-Recursive b(k):’ 
READ(*,100) FNAME 
OPEN(2,FILE=FNAME) 
READ(2,110) L 

write(*,*) Write k, b(k):’ 

DO 22 k=0,L 

READ(2,*) b(k) 

write(*,*) k,b(k) 

CONTINUE 

CLOSE (2) 

WRITE(*,*) ’Enter Filename for y[n]} Plot File:’ 
READ(*,100) FNAME 
OPEN(2,FILE=FNAME 
READ(2,100) TITL 
READ(2,110) NY 

READ(2,120) XQ 

READ(2,120) XQ 

DO 24 m=0,NY-1 

READ(2,120) y(m) 

CONTINUE 

CLOSE (2) 

WRITE(*,*) 

Setting up Plot File for z[n] 


WRITE(*,*) Enter Number of Points for z[n]: (le. 2048)’ 


READ(*,*) NZ 

WRITE(*,*) ’Enter time intervai (nanoseconds)’ 
READ(*,*) Tmax 

Tmin=0.0 

WRITE(*,*) Enter filename for z[n] response output’ 
READ(*,100) FNAME 

TITLE= ‘Inverse ARMA Filter z[n] Response’ 
OPEN(1,FILE=FNAME) 
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WRITE(1,100): TITLE 
WRITE(1,110) NZ 
WRITE(4,120):Tmin 
WRITE(1,120) Tmax 


i Cc Generating Inverse:Filter Coefficients. 
CG. Assuming b(0) .ne.-0 !!! 
Li=N 
‘ Ni=L 
DO 50 k=1,Ni 


50 ai(k)=-b(k)/b(0) 
bi(0) =1.0D00/b(0) 
DO 55 k=1,Li- : 
55 _ bi(k)=-a(k)/b(0) 
WRITE(*,*) 

C Initializing then iterating to form z[m] 
2(0)=bi(0)*y(0) 
write(1,120) z(0) 

DO 77 m=1,NZ-1 
z(m)=0.0 

Lmax=MIN(m,Li) 

DO 60 k=0,Lmax 
IF((m-k):GE:NY)-GO TO 60 
z(m)=z(m) +bi(k)*y(m-k) 

60 Continue 
Nmax=MIN(m,Ni) 
DO-66 k=1,Nmax 

‘i 66  z(m)=z(m)+<ai(k)*z(m-k) 

C — write(*,*) z(m) 

WRITE(1,120) 2(m) 
e 77. CONTINUE 
CLOSE(1) 





CALCULATING LATE-TIME ENERGY 


aQaa00 


WRITE(*,*) ENTER TL (Nsec) FOR LATE-TIME ENERGY CALCULATION’ 
WRITE(*,*) 
READ(*,*) TLATE 
WRITE(*,*) ENTER TE (Nsec) FOR LATE-TIME ENERGY CALCULATION’ 
WRITE(*,*) 
READ(*,*) TE 
NLATE=NINT((TLATE/Tmax)*NZ) +1 
NE=NINT((TE/Tmax)*NZ)-+1 
ENERGY=0.0 | 
DO 90 N=NLATE,NE 

90 ENERGY=ENERGY+Z(N)*Z(N) 





i ENERGY=ENERGY/(NE-NLATE +1) 
WRITE(*,*) 
WRITE(*,*) ‘LATE TIME ENERGY =’, ENERGY 
? WRITE(*,*) ‘Thats All Folks ! 
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400 
110 
120 


FORMAT(A) 
FORMAT(I5)- 
FORMAT(E12.6), 
STOP 

END 
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